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,  Existing  detailed  hourly  energy  analysis  programs  do  not 
adequately  model  the  heat  transfer  between  buildings 
and  the  ground.  A  simple  model  of  the  ground  heat 
transfer  compatible  with  both  existing  hourly  energy 
analysis  programs  and  simpler  building  models  is  vital 
as  energy  conservation  techniques  reduce  the  above¬ 
ground  heat  loss  and  building-ground  heat  transfer 
becomes  more  significant. 

This  study  extends  present  techniques  from  the  strictly 
geometric  context  of  the  numerical  solution  methods  to 
the  more  conceptual  environment  of  simplified  models. 
Specifically,  these  concepts  are  applied  to  the  problem 
of  heat  conduction  through  slab-on-grade  surfaces. 

Tested  over  a  broad  range  of  climatic  conditions,  the 
multiple-input  transfer  function  model  calculates  slab 
heat  flux.  The  accuracy  of  the  model  is  dependent 
upon  the  accuracy  of  the  input  data;  however,  some 
reasonable  approximations  to  the  necessary  input  data 
can  give  acceptable  results. 

The  full  capability  of  the  model  was  not  tested  in  this 
study.  Further  work  to  develop  a  definition  of  the  net¬ 
work  parameters  based  on  characteristic  length  could 
expand  the  use  of  the  model  to  nonsquare  and  possibly 
even  nonrectangular  surfaces.  .  ,  .1 
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MULTIPLE- INPUT  TRANSFER  FUNCTION  MODEL 
OF  HEAT  TRANSFER  FROM  SQUARE  SLAB  FLOORS 

1  INTRODUCTION 

Existing  detailed  hourly  energy  analysis  programs  such  as  BLAST’  do 
not  adequately  model  the  heat,  transfer  between  buildings  and  the 
ground.  Although  the  model  of  the  building  can  oe  very  complex,  the 
models  of  the  building-ground  heat  transfer  mechanisms  are  generally 
incongruously  simple.  BLAST,  for  example,  uses  a  one-dimensional 
response  factor  model  with  a  single  monthly  average  ground  tempera¬ 
ture  to  iefir.e  ail  building-ground  heat  transfer.  A  simple  model  of 
the  ground  heat  transfer  compatible  with  both  existing  hourly  energy 
analysis  programs  and  simpler  building  models  becomes  more  vital  as 
energy  censervat Lon  techniques  reduce  the  above-ground  heat  loss  and 
bu i Idi ng-grouna  heat  transfer  becomes  more  significant. 

Ceylan  and  Myers  [1|  developed  a  response -coefficient  method  for 
multidimensional  heat  conduction  problems  which  is  substantially  more 
efficient  than  finite-difference  or  finite-element  methods.  Addi¬ 
tionally,  it  provides  a  response  coefficient  model  of  the  system 
which  can  be  used  with  any  input  data  which  can  be  approximated  by  a 
continuous,  piecewise  linear  function.  Seem  [2]  developed  a  proce¬ 
dure  for  calculating  multidimensional  transfer  functions  which  elimi¬ 
nates  some  of  the  computationally  expensive  steps  of  the  Cev] an  and 
Myers  method.  These  multidimensional  methods  have  been  applied  to 
strict.]’/  geometric  heat  conduction  problems.  The  objective  of  this 
study  is  to  extend  these  techniques  from  the  strictly  geometric 
context  of  the  numerical  solution  methods  to  the  more  conceptual 
environment  of  simplified  models.  Specifically,  these  concepts  will 
be  applied  tc  the  problem  of  heat  conduction  through  slab-on-grade 
surface: . 

i  ’  ,  !  y  .  ye  <'■'  "i .  ■  ■  r  t-  ; -i  y  n  jp  i  cs  (UPAS’!  )  was  developed  by  'JSACERL  and 

■  •  "  :  iVi  fr.'.oi;  ! or  military  construction  projects. 
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2  CONCEPT 


Many  physical  systems,  including  thermodynamic  systems,  can 
be  approximated  using  lumped-system  analysis.  In  this 
approach  the  system  is  described  as  a  series  of  lumped,  lin¬ 
ear,  dynamic  elements  defined  by  ordinary  differential  equa¬ 
tions.  The  network  analogy  provides  a  simple  visualization 
of  this  concept.  In  a  network  model  of  a  thermal  system, 
temperatures  are  represented  by  nodes  with  a  lineal  tempera¬ 
ture  distribution  between  each  pair  of  nodes.  Physical 
properties  are  considered  to  be  uniform  between  each  pair  of 
nodes,  but  can  vary  from  pair  to  pair.  Energy  balance  equa¬ 
tions  are  written  for  each  node  and  the  system  of  equations 
solved  for  unknown  temperatures  and  heat  fluxes.  The 
validity  of  the  system  model  is  dependent  on  the  accuracy  of 
the  assumptions  of  uniform  temperature  at  each  node  and  lin¬ 
earity  between  nodes. 

Without  defining  specific  geometric  or  environmental  proper¬ 
ties,  the  matrices  forming  the  energy  balance  equations  of 
the  nodes  are  constructed  using  state  space  representation 
resulting  in  the  state  equation 

a  \  (i) 

_ =  q  \  +  25  £  ' 

dt 

and  the  output  equation 

o  =  cx  +  vr  (2) 

The  matrix  Y  contains  the  unknown  temperatures  (state 
variables).  I  is  the  matrix  of  known  temperatures  (input 
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variables) .  q  is  the  matrix  of  fluxes  (output  variables) . 
Matrices  ,  25,  t,  and  D  are  coefficient  matrices.  The 

size  of  the  matrices  and  the  values  of  the  elements  will  be 
determined  by  the  specific  model.  Once  the  coefficient 
matrices  are  defined  and  the  input  values  identified,  the 
first  order  differential  equations  can  be  solved.  The 
method  of  Seem  [2]  is  used  to  solve  the  system  of  equations. 
In  this  formulation,  the  time  series  of  input  variables  is 
modelled  as  a  continuous  piecewise  linear  function  by  the 
equation 


i  (0  =  r, 


(  T  -  0 
6 


x  *  6 


(3) 


Using  this  function  for  the  inputs,  the  differential 
equations  are  solved,  and  substituted  for  X  into  the 
equation  (2)  resulting  in  an  equation  relating  the  system 
outputs  to  the  system  inputs.  This  equation  is  known  as  the 
transfer  function  equation  and  is  of  the  form 


(■>) 


i-o 


i-o 


where 

Q,  =  vector  of  output  variables  (heat  flux)  at  time  i 

=  transfer  function  matrix  for  temperature  inputs 
at  time  j 

j  =  designator  identifying  a  point  in  time,  where  j=0  is  the 
current  time,  j=l  is  one  time  step  prior  to  the  current 
time  and  so  on 
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t  =  time  of  interest 
b  =  time  step 

(.!,  =  vector  of  input  variables  (known  temperatures) 
at  time  i 

e  =  scalar  constant  for  adjusting  the  effect  of  previous 
outputs  on  the  output  at  the  time  of  interest. 


A  transfer  function  is  defined  as  the  ratio  of  the  output 
variables  of  a  system  in  state  space  to  its  input  variables 
(also  in  state  space) .  In  this  way,  the  transfer  function 
represents  the  dynamics  of  a  linear  time-invariant  system. 
The  transfer  function  matrices  are  dependent  on  the  system 
and  inputs,  but  only  on  the  functional  form  of  the  inputs; 
therefore,  any  input  which  can  be  adequately  modelled  by  the 
continuous  piecewise  linear  function  noted  above  can  be  used 
with  the  transfer  function  matrices  to  model  its  effect  on 
the  system  This  is  particularly  useful  in  the  modelling  of 
building  systems  where  the  input  conditions  of  climate  are 
unpredictable  and  highly  variant  in  time  and  geographic 
location. 
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3  MODEL  DEVELOPMENT 


edge  region  temperature  (/J  ,  the  daily  average  ground  sur¬ 
face  temperature  (/'/)  and  the  deep  ground  temperature  (T  d)  . 
The  3  state  variables,  the  temperatures  at  the  remaining 
nodes  (7',,  /  2,  /  3)  ,  are  allowed  to  float  and  consequently 

have  some  thermal  capacitance  attributed  to  them.  The  tem¬ 
perature  nodes  are  related  to  each  other  as  shown  in  the 
figure.  Between  attached  pairs  of  temperature  nodes,  there 
exists  some  thermal  resistance.  The  definition  of  these 
resistances  and  capacitances  is  discussed  in  Section  3.2. 


3.2  BASIC  EQUATIONS 

BASIC  EQUATIONS  Energy  balance  equations  are  written  for 
each  node  resulting  in  4  state  equations  of  the  form 


dl\  f  (5) 


for  i  =  1  to  4,  and  3  state  output  equations  of  the  form 


/ 

/-  l 


for  i  -  1  to  3 ,  where 
C ,  =  thermal  capacitance  at  node  V, 

G, ,  = —  =  inverse  of  the  thermal  resistance  between 

h  n 

nodes  I  ,  and  V, 

These  7  equations  can  be  written  more  conveniently  in  matrix 
form 


6 


(7) 


d_\ 

dl 


AX +  31/ 


and 


Q  =  C  X  +  'DU 


(8) 


where 
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The  coefficient  matrices  A,  B ,  C,  and  ©  define  the 
relationships  of  all  temperature  regions  in  the  system  to 
all  others.  They  involve  geometric  factors  such  as  the  area 
through  which  heat  is  transferred  from  one  region  to 
another,  and  physical  properties  such  as  the  density  and 
thermal  conductivity  of  various  regions.  The  goal  of 
defining  the  elements  of  the  coefficient  matrices  is  to  make 
it  possible  to  generate  transfer  function  equations  for  any 
system  from  its  basic  physical  parameters  rather  than  as  is 
frequently  done  in  electro-mechanical  systems  -  by  testing 
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the  system  itself.  Because  the  important  aspect  of  the 
equations  is  the  thermal  relationships  between  regions,  the 
model  is  not  strictly  geometric. 


The  first  step  in  defining  the  matrix  coefficients  is  iden¬ 
tifying  the  properties  which  make  up  the  elements  of  C  and 
C.  The  basic  form  allows  for  the  description  of  several 
heat  transfer  mechanisms  given  the  appropriate  temperatures. 
For  conduction,  the  equation  becomes 


dT 

Q  =  kA  — 
d  x 


(17) 


or,  in  the  spatially  discretized  form  used  for  this  model 


(18) 


kA 


In  this  case  Clf  is  defined  as  the  conductance,  where 


kK,  =  thermal  conductivity  applicable  to  the  volume 
between  nodes  i  and  j 

/!,,  =  cross-sectional  area  through  which  heat  is  trans¬ 
ferred  between  nodes  i  and  j 

I-,,  =  distance  between  nodes  i  and  j. 

Although  this  model  does  not  contain  convective  or  radiative 
heat  transfer,  these  mechanisms  can  be  supported  by  the 
model  by  setting 


9 


G  =  h  A 


(19) 


for  convective  heat  transfer,  and 


C  =  lirA 


for  radiative  heat  transfer,  where 


h  =  convective  heat  transfer  coefficient 


hr  =  the  effective  linearized  radiative  heat  transfer 
coefficient. 

The  thermal  capacitance  C  is  derived  from  the  transient 
equation 


,  dT 
<3-pc„i  -j7 


so  that 


Ci  =  P.Cpjl-'i 


where 


p,  =  density  of  the  region  of  soil  at  T, 

(' ,,,  =  specific  heat  of  the  region  of  soil  at  7, 
l  ,  =  volume  of  the  region  of  soil  at  Tt. 
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Both  thermal  conductance,  GtJ,  and  thermal  capacitance,  Ct> 

are  composed  of  geometric  factors  ( Lit ,  Au,  and  V ()  as  well 
as  soil  properties  (ku,  p, ,  and  cp.)  .  These  will  be  dis¬ 
cussed  separately. 


3.3  GEOMETRY 


Bahnfleth's  study  of  undisturbed  ground  temperature  patterns 
[3]  shows  two  distinctly  different  zones  of  temperature 
fluctuation:  a  relatively  fast  zone  near  the  ground  surface 
where  the  temperature  changes  are  in  scale  with  the  tempera¬ 
ture  changes  of  the  forcing  temperature,  and  a  slower  zone 
where  temperature  fluctuations  are  strongly  damped.  Because 
the  response  rate  of  the  near-surface  zone  is  more  similar 
to  the  response  rate  of  typical  building  components  than  to 
that  of  the  remainder  of  the  earth,  it  was  decided  to  model 
the  near-surface  earth  and  the  remainder  of  the  earth  as 
attached  but  distinct  components.  The  point  of  separation 
of  these  zones  is  the  diurnal  penetration  depth,  or  roughly 
0.5  meters  below  the  surface.  The  temperature  at  this  point 
remains  nearly  constant  over  a  day  at  the  daily  average 
ground  surface  temperature. 

Horizontal  maps  of  ground  temperature  beneath  buildings  show 
a  circular  pattern,  conseguently ,  a  cylindrical  coordinate 
system  was  used  to  produce  an  axisymmetric  two-dimensional 
model.  Horizontal  temperature  nodes  are  set  at  the  slab 
center,  the  edge-equivalent  radius,  and  the  location  where 
the  ground  temperature  is  unaffected  by  the  building  (far- 
field) .  The  edge-equivalent  radius  is  calculated  as  the 
radius  of  a  circle  of  equivalent  slab  perimeter,  or 
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(23) 


P 

'  11  ~2r 

In  this  fashion,  a  10m  by  10m  slab  is  mapped  to  a  circle  of 
radius 


r 


p 


12 
2  u 


6.37  m 


(24) 


Vertical  temperature  nodes  are  set  at  the  diurnal  penetra¬ 
tion  depth  of  the  surface  temperature  wave  (approximately 
u.5m  below  the  surface),  the  annual  penetration  depth 
(approximately  15m  below  the  surface)  and  the  depth  of  the 
point  of  inflection  or  knee  of  the  undisturbed  temperature 
profile  (approximately  4m  below  the  surface).  Studies  ([3], 
[4],  [5])  of  underground  temperature  patterns  show  a  shape 

which  could  be  approximated  by  linear  temperature  profiles 
between  these  temperature  nodes  (Figure  2) . 

The  area-eguivalent  radius  is  used  for  calculations  in  the 
vertical  plane.  It  is  calculated  as  the  radius  of  the 
circle  which  has  the  same  area  as  the  slab,  i.e. 


n  (25) 

r,'  =  \/  Ti  • 

Therefore,  for  calculations  in  the  vertical  direction,  the 
10m  by  10m  slab  is  mapped  to  a  circle  of  radius 


/  100 

- =  3.6  I m 

V  n 


(26) 


12 


so  that  area  is  preserved. 


Figure  2:  Undisturbed  Ground  Temperature  Profile 


D 

e 

P 

t 

h 
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Although  this  model  cannot  be  reproduced  graphically,  it 
accounts  for  both  the  perimeter  and  area  effects  of  ground- 
coupled  heat  transfer. 

These  geometric  relationships  were  used  to  build  the  geomet¬ 
ric  matrices  which  were  used,  in  turn,  to  develop  the  coef¬ 
ficient  matrices. 
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Three  matrices  constitute  the  geometric  factors  of  the 
1,  L,  and  V.  Matrix  A  is  the  symmetric  7x7  matrix 


A  l? 

-1,3 

/1  1  d 

/iiA 

f  '1,2 

•  ??. 

'1  23 

■4*b 

12/ 

"'1  2d 

•■1,3 

''1  23 

•  1  33 

'1  3  b 

'*1  3  / 

71  3  d 

'i3e 

•'1  1  b 

2  b 

'  1  3  b 

''1  bb 

1  «>/ 

1  bd 

'‘1  bt> 

I  .  / 

/l  2f 

>3/ 

'V 

■'«// 

'I 

\ 

''1  ?i  1 

•'1  3d 

1  ()</ 

1  /<! 

1  dd 

•'Id, 

i 

\  lln 

1  2p 

■1.,, 

1  dp 

w 

where  the  subscripts  refer  to  the  path,  that  is,  At,  is  the 
area  through  which  heat  is  transferred  from  /',  to  T  r 


Matrix  /  is  a  similar  symmetric  7x7  matrix: 


flu 

/,2 

In 

/•  1  b 

/,/ 

/  1  d 

!  v; 

/  2?. 

23 

/•2  b 

1-  2  f 

/-  2d 

l  ,3 

I  23 

/•  33 

/•  Jb 

/  3  / 

/•3  d 

£3,  ! 

/•lb 

/•2b 

/■  3b 

/•  bb 

l-bf 

/'  bd 

be 

/  .  / 

1-2.1 

/'  3  / 

/  b/ 

hi 

/•  fd 

^  | 

1  2d 

/  3d 

/•  bd 

/•  /d 

/•  dd 

/do  i 

W 

1-2, 

/■J, 

/  bp 

/'  dp 

For  the  7  node  network  model,  not  all  of  the  nodes  are 
connected,  therefore  the  elements  which  relate  unconnected 
nodes  to  each  other  become  0  leaving 
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/  o  /ll2  0  /llb  0  Ald  o  \  (29) 

12  0  A?:i  0  ^  2d  ^2e 

0  A  23  0  0  /l3y  ^3d  0 

,4  =  /1,k  0  0  0  0  0  Ab, 

0  0  /13,  0  0  0  A,e 

,  A  \  d  A  2(i  A  3cl  0  0  0  0  I 

\  0  a2,  o  ,i„.  /I /e  0  0  / 

and 

/()  /.,  2  0  Z.lb  0  Z.ld  0\  (30) 

(-12  0  Z.  ?3  0  0  Z.  2d  (■  2e 

0  (-23  0  0  Z  3y  L  3d  ® 

/.=  Zlb  0  0  0  0  0  zbe  . 

0  0  Z3/  0  0  0  Z/e 

Z-,d  (-2d  (-3d  0  0  0  0 

\  0  (-2,  0  (-6,  (-/„  0  0  / 

l  is  the  vector 


The  definition  of  the  elements,  or  network  parameters, 

1(/,  Z.t/,  and  l',  is,  in  part,  independent  of  the  model 

structure,  but  is  based  on  the  geometry  presented  above  and 
is  described  further  in  Section  4. 
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3.4  SOIL  PROPERTIES 


Soil  properties  are  represented  by  another  7x7  symmetric 
matrix: 


/  0  fr.2 

/  A|,  0 

'  0  krj 


0  0 


VA'l,f  A'  ?./ 
0  A-., 


0  A,h  0 

A'  7;,  0  0 

0  0  A  t/ 

0  0  0 

A'  t/  0  0 

A  0  0 

o  A,„.  A-/(, 


A-,d  0  \ 

k  ;  p  \ 

o 

o  A  ,p 

0  A'/c 

0  0  I 

0  0  / 


and  two  vectors, 


P  =  ( P  :  P . ) 


and 


(32) 


(33) 


)• 


(34) 


Thermal  properties  of  the  soil  can  be  defined  separately  for 
each  energy  balance  equation.  Individually,  these  equations 
assume  constant  thermal  properties.  Consequently,  the 
properties  are  defined  as  the  ’'effective"  value  of  the 
thermal  properties  in  the  region  specified  by  the  equation. 
The  need  for  an  "effective"  conductivity  is  described  more 
thoroughly  in  Bahnfleth  [3].  Effective  conductivity  is 
defined  as 


A 


P 


(35) 
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where 


lk  =  the  thickness  of  layer  k,  and 
hk  =  the  conductivity  of  layer  k. 

Calculated  from  this  equation,  kt/  is  the  effective  value  of 

the  soil  conductivity  in  the  region  through  which  heat 
transfers  between  T,  and  T 


3.5  INPUTS 


The  inputs  to  the  transfer  function  equation  are  (referring 
to  Figure  1)  the  temperatures  of  the  slab  core  area  (T b)  and 
the  slab  edge  area  (Te)  ,  and  the  undisturbed  ground  tempera¬ 
tures  near  the  surface  (7  /)  and  in  the  deep  ground  {T d)  . 
Because  the  "top"  nodes  are  defined  at  the  diurnal 
penetration  depth,  their  temperatures  can  be  approximated  by 
the  daily  average  of  the  surface  temperature,  i.e.  Tb  =  the 
daily  average  slab  center  temperature  and  T ,  =  the  daily 

average  ground  surface  temperature.  The  undisturbed  ground 
surface  and  deep  ground  temperatures  can  be  determined  a 
priori  using  a  variety  of  algorithms. 
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4  NETWORK  PARAMETER  SPECIFICATION 


4.1  METHOD 

The  network  parameters  are  the  elements  which  compose  the 
geometric  matrices  /I,  / ,  and  I  .  Their  definition  depends 

on  the  method  of  discretizing  the  geometry  of  the  system.  In 
other  words,  the  magnitude  of  the  element  is  dependent  on 
the  sizes  of  the  regions  assumed  to  be  at  the  specified 
temperatures.  Because  few  of  the  regions  are  actually  iso¬ 
thermal,  the  allocation  of  area  and  volume  to  a  specific 
temperature  must  be  based  on  some  method  or  algorithm. 

In  order  to  evaluate  the  validity  of  the  postulated  network 
parameters  a  test  system  which  fixes  input  conditions  and 
environmental  parameters  was  established.  A  base  case  using 
the  same  fixed  conditions  was  used  to  compare  both  the  total 
flux  data  and  the  form  of  the  daily  average  flux  curve. 

This  provided  not  only  data  on  the  accuracy  of  the  model  but 
also  clues  to  the  nature  of  any  inaccuracies  that  might 
occur. 


4.2  BASE  CASE 

The  baseline  for  evaluation  of  the  model  is  a  detailed 
finite  difference  model  ( FDM)  of  heat  transfer  from  slab 
floors  developed  by  Bahnfleth  [3].  A  three-dimensional 
model  of  a  slab-on-grade  and  the  soil  beneath  it  is  solved 
by  numerical  techniques  in  the  program  SLAB3D.  The  space 
above  the  slab  is  defined  by  a  constant  room  air  tempera¬ 
ture.  Undisturbed  soil  temperature  distributions  as  calcu¬ 
lated  from  subroutine  TEARTH  are  used  as  the  far-field 
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boundary  temperatures.  The  deep  ground  boundary  can  be 
either  a  specified  flux  or  a  specified  temperature  plane. 
Bahnfleth  generated  data  for  a  variety  of  surface  shapes  and 
sizes  as  well  as  a  number  of  diverse  climatic  conditions. 


4.3  DESCRIPTION  OF  TEST  SYSTEM 

4.3.1  GEOMETRY 

Because  the  model  was  developed  based  on  a  square  slab-on- 
grade,  the  development  of  the  network  parameters  was  based 
on  this  geometry.  Base  case  data  were  available  for  two 
sizes  of  square  slabs,  a  12m  x  12m  square  and  a  45m  x  45m 
square.  The  flux  per  unit  area  calculated  by  the  finite 
difference  model  for  these  two  slabs  for  a  calendar  year  in 
Minneapolis  MN  is  shown  in  Figure  3.  Due  to  the  dominance 
of  the  edge  effect  in  the  smaller  slab,  the  annual  flux 
variation  is  much  greater  for  the  smaller  slab.  If  the 
effect  of  the  balance  between  the  perimeter  loss  and  edge 
loss  is  to  be  accommodated,  it  is  important  that  the  test 
system  strongly  exhibits  this  effect.  Consequently,  the 
smaller  12m  x  12m  slab  was  used  as  the  primary  cest  system. 


Figure  3:  FDM  of  Two  Square  Slabs 
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4.3.2  SOIL  PROPERTIES 

Although  the  system  model  can  support  variable  soil  proper¬ 
ties,  this  capability  is  not  tested  in  this  study.  Avail¬ 
able  base  case  heat  flux  data  were  calculated  using  constant 
soil  properties.  The  same  properties  are  used  in  this 
study.  They  are: 

p  =  1  200  kg /in3 


c  p  =  1 200  J/kg/K 

k  =1.00  U'/m/k 

4.3.3  INPUTS 


For  this  study,  existing  data  from  Bahnfleth's  one¬ 
dimensional  semi-infinite  solid  model  of  the  heat  transfer 
in  undisturbed  earth  [3]  were  used  to  provide  input  data  for 
the  nodes  at  the  far-field  and  the  deep  ground.  The  hourly 
ground  surface  temperatures  calculated  by  this  model  vary 
with  local  climatic  conditions  while  the  deep  ground  temper¬ 
ature  is  constant  at  the  annual  average  ground  surface  tem¬ 
perature.  Because  the  far-field  node  is  placed  at  the 
diurnal  penetration  depth,  the  daily  average  ground  surface 
temperature  rather  than  the  hourly  ground  surface  tempera¬ 
ture  is  used  as  input  to  the  multiple  input  transfer  func¬ 
tion  model. 
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Exact  data  for  the  daily  average  slab  center  and  slab  edge 
temperatures  were  not  available  from  the  base  case  data  set. 
This  is  the  usual  case  with  most  energy  analysis  programs. 
Therefore,  although  it  is  possible  in  the  network  mode  to 
include  a  temperature  difference  between  the  slab  center  and 
the  slab  edge,  daily  average  slab  surface  temperatures 
(which  are  equal  to  the  ground  temperature  at  the  diurnal 
penetration  depth)  were  used  for  both  these  temperatures. 

The  network  parameters  developed  using  the  assumption  of  an 
isothermal  slab  should  be  appropriate  for  use  with  energy 
analysis  programs,  such  as  BLAST  and  DOE-2  which  use  the 
same  assumption. 

Data  were  available  from  the  base  case  data  set  for  four 
locations,  Minneapolis  MN,  Medford  OR,  Philadelphia  PA,  and 
Phoenix  AZ.  These  data  were  generated  originally  from  Typi¬ 
cal  Meteorological  Year  (TMY)  weather  data  for  these  sites. 
Figures  4,  5,  6,  7  show  the  daily  average  air  temperatures 
at  these  four  locations. 

In  order  to  develop  the  most  responsive  model  network  param¬ 
eters,  the  most  rigorous  weather  conditions  were  used.  The 
figures  indicate  that  the  weather  data  for  Minneapolis  MN 
would  provide  the  most  demanding  conditions  for  the  model. 

In  addition  to  a  large  annual  temperature  variation,  the 
temperature  variation  from  day  to  day  is  also  larger  in  the 
Minneapolis  data  than  in  that  of  the  other  three  locations. 
During  the  development  of  the  network  parameters,  therefore, 
the  data  derived  using  Minneapolis  weather  were  utilized  to 
evaluate  the  accuracy  of  the  model  fit.  The  model  was  then 
tested  later  at  the  remaining  three  locations. 
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Figure  6:  Daily  Average  Air  Temperature  —  Philadelphia  PA 


Figure  7:  Daily  Average  Air  Temperature  —  Phoenix  AZ 


4.4  GEOMETRIC  DEFINITION  OF  NETWORK  PARAMETERS 


The  geometric  definition  of  the  network  parameters  is  based 
on  a  series  of  assumptions  about  the  size  and  shape  of  the 
temperature  regions  beneath  the  slab  as  well  as  the  geometry 
of  the  slab  itself.  The  geometry  of  the  slab  and  the  ground 
in  the  vicinity  of  the  slab  is  the  foundation  for  these 
assumptions.  For  the  test  system,  the  slab  geometry  is 
defined  by 


A  =  I2.v  12  =  II 4  in 


l\  =  4 a  12  =  18m 


r-"VT"6*/7m 

rs 

r  =  —  =  7.64  m. 

'  2  k 


(36) 


The  ground  geometry  has  already  been  fixed  by  the  network 
model  [see  Figure  1]. 

The  elements  of  the  matrix  /.  can  then  be  assigned  as  the 
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Referring  to  Figure  1,  note  that 


L  lb  ~  !-2e  ~  D  x 

L  I  d  =  ^2d  =  l  3d  =  ®2 


L  fe  1-23  D3 


and 


'  be 


~  L  12  1  P  ■ 


Substituting  Dt.  D?,  D 3,  and  rp  into  L, 
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0i 

0 

l  D  2 

\  0 


rp  0  D, 
0  D3  0 
D3  0  0 
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0  9,0 
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The  elements  of  /I  are  defined  as  the  area  through 
heat  is  transferred  between  /  ,  and 
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(39) 

(40) 
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(43) 
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The  determination  of  the  elements  of  A  is  more  complex  than 
the  determination  of  the  elements  of  L.  Although  the  model 
is  not  strictly  geometric,  the  heat  transfer  areas  can  be 
initially  postulated  based  on  geometric  considerations. 
First,  it  is  assumed  that  Alb  =  Ald,  A2e  =  A2d,  and  A3/  =  A3d. 
Then  the  area  through  which  heat  is  transferred  between 
T b  and  T ,  can  be  estimated  by  the  area  of  the  slab  which 
can  be  approximated  by  the  daily  average  slab  surface  center 
temperature,  Tb.  Studies  [3], [4],  and  [5]  have  shown  that 
temperature  gradients  across  horizontal  ground-contact  sur¬ 
faces  are  small  over  most  of  the  surface  and  relatively 
large  near  the  edge.  For  this  model  the  area  in  the  plane 
of  the  ground  surface  considered  to  be  at  Tb  is  the  area  of 
the  slab  minus  the  area  near  the  slab  edge,  or 

nr.-d,)*  <44> 

where  d,  is  the  distance  from  the  slab  edge  to  the  location 
where  the  slab  surface  temperature  is  approximately  T b . 

The  area  through  which  heat  is  transferred  between 
Te  and  7 2  can  be  defined  as  the  area  where  the  edge  effect 
predominates.  This  area  is  specified  as  the  area  within  a 
definable  distance,  d2  of  the  building  edge  minus  the  area 
already  associated  with  7b,  or 

l2.=  'l2rf-w(ra  +  rf2)2-.-l16  (45> 

where  d-.  is  the  distance  the  edge  effect  region  extends 
beyond  the  slab  edge. 
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The  remaining  area  is  associated  with  the  far-field  tempera¬ 
ture,  T It  is  calculated  from  the  equation, 

A  3f  ~  A3d  =  Tl(r  a  +  D ~Au>~A2a  (46) 

where  D3  is  the  distance  from  the  edge  of  the  slab  to  the 
undisturbed  ground  in  the  far-field,  or  12.5m  [3], 

The  area  through  which  heat  is  transferred  between 
7",  and  V ,  is  calculated  from  the  Fourier  equation  for  con¬ 
duction  through  a  hollow  cylinder 

f 2n/A  (47) 


By  the  convention  of  the  matrix  definition,  all  equations 
are  cast  in  the  general  form 
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Therefore, 
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so  that 
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(50) 
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2  It/.  2e^-e/  2n.(D])(D3) 


Similarly, 
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(51) 


(52) 


It  is  postulated  that  the  slab  is  isothermal  to  within  c/ , 

the  slab  edge.  It  is  therefore  estimated  that  the  area 
through  which  heat  is  transferred  between  T b  and  T„  can 
also  be  calculated  correspondingly: 


lb,= 


2ii/.|t,(c/,)  _2n(/J,)(d1) 
1  n  r— ~  In 


(53) 


Again,  Al2  is  calculated  in  a  similar  fashion: 
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(54) 


Of 
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The  volumes  associated  with  temperatures  T,t  T2,and  T 3  are 
calculated  by 

i' i i  <55> 

1^2 

and 

l  3  =  ^  3/  J 

I)  i ,  D2.  and  /7;)  are  defined  by  temperature  profiles  of  the 

undisturbed  ground  and  do  not  change  with  slab  size  or 
shape.  Referencing  Section  3.3  and  Figure  2,  D,  is  the 
distance  between  the  inflection  point  and  the  diurnal 
penetration  depth  or, 

D,  =  4.0m  -  0.5m  =  3.5m.  (56) 

D 7  is  the  distance  between  the  undisturbed  deep  ground  and 
the  inflection  point  or, 

D?  =  15m  -  4  m  =  11m  (57) 

and  U3  is  the  distance  from  the  edge  of  the  slab  to 
undisturbed  ground  temperature  in  the  vertical  plane, 

D3  =  12.5m.  (58) 

h  |,  h2.  and  h3  are  assumed  to  be  equal  to  each  other  and 
equal  to  half  of  the  depth  of  the  entire  system: 

/i,  =  /i2  =  /ia  =  14.5/2  =  7.25m.  (59) 
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The  results  of  Bahnfleth  [3]  show  that  for  rectangular 
slabs,  surface  temperature  is  nearly  constant  to  within  1.0m 
of  the  edge  regardless  of  slab  size  or  aspect  ratio.  Based 
on  these  data  d,  is  assigned  the  value  of  1.0m. 

In  the  edge  effect  region,  the  direction  of  heat  flux 
changes  dramatically  over  the  year.  Temperature  profiles 
developed  by  Bahnfleth  [3],  and  Kusuda  [5]  show  that  in  the 
region  between  the  slab  edge  and  roughly  3m  outside  the 
building  edge,  the  direction  of  heat  flux  varies  dramati¬ 
cally  over  the  annual  cycle.  It  is  desirable  that  this 
region  be  distinguished  from  the  regions  where  the  heat  flux 
patterns  are  more  consistent  over  the  annual  cycle.  There¬ 
fore,  this  region  is  characterized  as  the  edge  region  and 
its  area  (in  the  plane  of  the  ground  surface),  A?e,  is 
encircled  by  a  boundary  approximately  3.0m  beyond  the  edge 
of  the  slab.  Therefore,  d2  is  assigned  the  value  of  3.0m. 

For  the  test  system,  /.  becomes 
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Substituting 

/' 

ra. 

<7  i 

and  d 

2  into  the 

area  equations 

and  introducing  the  values  into  the  matrix  A  gives 
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The  volume  vector  then 


becomes 


v  = 


(62) 


4.4.1  MULTIPLE-INPUT  TRANSFER  FUNCTION  COEFFICIENTS 


Program  GTF  [see  Appendix  C]  calculates  the  multiple-input 
ground  transfer  function  (GTF)  coefficients  and  scalar  con 
stants  from  the  model  structure  and  geometry  and  soil  prop 
erties  using  the  method  of  Seem.  The  GTF  coefficients  and 
scalar  constants  calculated  for  the  test  system  with  the 
geometric  network  parameters  are: 


S0 
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s ,  * 
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p  ,  =  -  2.9792  I  (64) 

p3  =  -0.9/9329. 


4.4.2  INITIAL  HEAT  FLUX  CALCULATIONS  USING  MULTIPLE  INPUT 
TRANSFER  FUNCTIONS 

Program  QCALC  (see  Appendix  C)  uses  the  GTFs  and  the  input 
temperatures  with  Equation  (5)  to  calculate  the  daily  aver¬ 
age  heat  flux.  The  results  are  divided  by  the  slab  area  in 
order  to  be  compatible  in  units  with  the  base  case  data. 
Figure  8  shows  the  flux  plot  of  the  GTF  model  vs  the  finite 
difference  model.  Figure  9  plots  the  daily  average  flux  of 
the  GTF  model  and  its  difference  from  the  FDM .  Although 
annual  period  of  the  flux  curve  appears  nearly  correct,  its 
amplitude  is  much  too  small.  Because  the  magnitude  of  the 
maximum  flux  is  fairly  accurate,  it  is  likely  that  the  pri¬ 
mary  cause  of  the  error  is  the  assumption  of  too  much  mass 
in  the  system.  The  network  parameters  were  based  on  the 
assumption  that  the  mass  associated  with  the  system  is  com¬ 
posed  of  three  cylinders  of  equal  depth  with  cross-sectional 
areas  equal  to  the  area  through  which  heat  is  transferred 
vertically  in  their  respective  regions.  The  resulting  vol¬ 
ume  vector  (Equation  (62))  heavily  weights  the  system  mass 
with  the  far-field  undisturbed  ground  temperature.  The 
result  is  that  the  mass  of  soil  directly  beneath  the  slab 
has  the  least  effect.  This  is  obviously  incorrect.  Various 
changes  to  the  model  were  introduced  to  improve  its  behav¬ 
ior.  These  are  discussed  in  Section  4.5. 
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DAILY  AVERACE  HEAT  FLUX  (W/mA2)  DAILY  AVERAGE  HEAT  FLUX  [W/m*2] 


Figure  8 :  Flux  —  FDM  and  Run  1A 
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Figure  9:  Flux  and  Difference  —  Run  1A 
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4.5  PARAMETER  REFINEMENT  USING  EMPIRICAL  METHODS 


4.5.1  MODEL  BASED  ON  12m  X  12m  SLAB 

The  effect  of  soil  volume  on  the  shape  of  the  annual  flux 
curve,  although  evident,  is  not  easily  described  by  geomet¬ 
ric  techniques.  It  may  be  possible  to  improve  the  accuracy 
of  the  model,  particularly  in  cases  where  the  daily  average 
flux  per  unit  area  is  small,  by  adjusting  the  model  parame¬ 
ters  based  on  improved  fit  to  the  base  case  data. 

The  original  parameter  set  based  on  system  geometry  showed 
an  excess  of  thermal  mass  in  the  system.  Therefore,  the 
volume  vector  was  revised  by  reducing  the  depth  of  the  cyl¬ 
inders  of  mass  associated  with  1'2  and  l’:i  to  lm  from  7.25m 
while  leaving  the  mass  associated  with  I  t  intact.  The 
volume  vector  then  became 


(7  58  \ 
l'  =  19b 

\86 7  ) 


in 


3. 


(65) 


The  remaining  parameters  were  held  the  same  as  the  original 
run.  Table  1  Run  IB  presents  the  parameter  set. 

Program  GTF  was  rerun  with  this  new  set  of  parameters  and 
the  results  are  shown  in  Figures  10  and  11. 
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Table  1:  Parameter  Sets  for  Runs  1A  -  1C 


RUN  1A 

RUN  IB 

RUN  1C 

AREA  [m2] 

144 

144 

144 

PERIMETER  [m] 

48 

48 

48 

A/P  [m] 

3.0 

3.0 

3.0 

157 

157 

157 
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493 

493 

^  be 
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284 

284 

'‘I  23 

891 

891 

891 

A  1  6  =  1  d 

105 

105 

105 

■T  —  w  - - 

195 

195 

195 

^3/  ^3d 

867 

867 

867 

758 

758 

758 

^2 
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I7  3 

6284 

867 
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Figure  10:  Flux 


FDM  and  Run  IB 


Figure  11:  Flux  and  Difference 


Run  IB 


Figure  10  plots  the  daily  average  heat  flux  over  the  cycle 
for  the  base  case  FDM  and  the  GTF  Run  IB.  Comparison  with 
Figure  8  shows  the  substantial  improvement  in  fit.  Figure 
11  shows  the  flux  and  the  difference  from  the  FDM  results  of 
Run  IB.  Comparing  Figure  9  with  Figure  11  it  can  be  seen 
that  model  IB  vastly  reduces  the  model  error  in  the  summer¬ 
time  data. 


The  new  set  of  parameters  leads  to  a  better  match  of  the 
base  case  data  from  which  it  can  be  inferred  that  the  pri¬ 
mary  flaw  in  the  original  parameter  set  was  indeed  an  excess 
of  mass.  Even  though  the  revised  set  of  parameters  is 
substantially  better  than  the  original  set,  the  figures 
indicate  that  more  improvement  is  required.  The  annual 
amplitude  of  the  curve  is  still  slightly  too  small  and  the 
new  set  of  parameters  presents  a  new  flaw  in  the  model:  the' 
flux  calculated  for  the  end  of  the  year  becomes  increasingly 
too  small.  This  seems  to  indicate  an  imbalance  in  the  mass 
distribution.  The  first  step  towards  resolving  this  situa¬ 
tion  was  to  reduce  the  mass  associated  with  the  far-field 
region  in  an  effort  to  both  reduce  the  total  mass  of  the 
system  and  place  more  emphasis  on  the  mass  immediately 
beneath  the  slab.  The  depth  of  the  cylinder  of  mass  asso¬ 
ciated  with  the  far-field  (fr3)  was  reduced  from  lm  to  0.1  m. 
resulting  in  the  new  volume  vector 


V  = 


f  758  \ 
195 

V  87  ) 


m 


(66) 


Once  again,  the  remaining  network  parameters  were  held  the 
same  as  the  original  run  (see  Table  l  Run  1C) . 
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The  new  GTFs  and  scalar  constants  were  used  to  recalculate 
the  daily  average  heat  flux.  The  resulting  flux  and  differ- 
e.'.ce  plots  are  given  in  Figures  12  and  13.  On  a  plot  of  the 
daily  average  flux  (Figure  12)  Run  1C  and  the  FDM  are 
practically  indistinguishable.  Figure  13  presents  the  flux 
calculated  in  Run  1C  and  the  difference  between  Run  1C  and 
the  FDM.  Except  for  a  few  points  early  in  the  annual  cycle 
when  the  flux  is  high,  the  difference  between  Run  1C  and  the 
FDM  is  less  than  1  W/nT2. 

The  match  between  the  FDM  and  Run  1C  is  significantly  better 
than  both  Run  la  and  Run  lb.  Table  2  gives  a  numerical 
comparison  of  the  three  GTF  models  with  the  base  case  finite 
difference  model.  It  can  be  seen  that  90%  of  the  data  were 
within  15%  of  the  FDM  and  the  percent  error  in  total  energy 
consumption  over  the  annual  cycle  is  -1.0%. 


Figure  12:  Flux  —  FDM  and  Run  1C 
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4.5.2  MODEL  CORRECTIONS  BASED  ON  45m  X  45m  SLAB 


Using  the  same  parameter  equations  as  Run  1C,  changing  only 
r p  and  r„,  the  GTFs  and  scalar  constants  were  recalculated 
for  the  larger  (45m  x  45m)  slab.  The  resulting  parameter 
set  is  shown  in  the  column  labeled  Run  2A  in  Table  3. 

Table  3:  Parameter  Sets  for  Runs  2A  -  2D 


RUN  2  A 

RUN  2  B 

RUN  2C 

RUN  2D 

AREA  l  rn  2  ] 

2025 

2025 

2025 

2025 

PERIMETER  [m] 
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180 

180 

180 
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11.25 

11.25 

11.25 
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Figure  14  compares  the  output  for  the  45m  x  45m  slab  of  the 
GTF  model  and  the  FDM  base  case. 
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Figure  14:  Flux  and  Difference  —  Run  2A 


There  is  a  distinct  droop  of  the  GTF  during  the  winter  sea¬ 
son,  the  time  when  the  edge  effect  is  strongest.  During  the 
summer  season  when  the  edge  effect  is  small,  the  GTF  model 
is  quite  accurate.  This  suggests  that  for  larger  slabs  the 
method  for  calculating  the  parameter  set  does  not  account 
adequately  for  the  edge  effect.  Returning  to  the  original 
energy  balance  on  the  edge  node 
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(67) 


and  recalling  that  for  this  study  Ta  =  Tbl  we  see  that  there 

are  three  likely  causes  for  underprediction  of  e^ge  loss. 
The  first  possibility  is  that  AZa  is  too  small  so  that  loss 
to  the  ground  beneath  the  edge  is  underpredicted.  An 
adjustment  to  A?a  should  produce  a  generally  linear  change 
in  the  flux,  although  the  seasonal  change  of  T 2  will  modify 
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that  effect  to  some  extent.  The  second  possibility  is  that 
the  area  A  fa  is  too  small  so  that  heat  flux  directly  to  the 
ground  surface  is  insufficient.  Thirdly,  the  temperature 
node  T 2  may  be  too  high.  This  is  most  likely  to  result  from 
a  deficient  A22  which  would  lead  to  a  reduction  in  the  heat 
transferred  from  the  node  beneath  the  slab  edge  (node  2)  to 
the  node  at  4  meters  depth  in  the  far-field  (node  3)  thereby 
maintaining  an  inappropriately  elevated  T 2. 

Two  studies  were  performed  in  order  to  determine  the  most 
likely  cause  of  the  underprediction  of  edge  loss  in  the 
larger  slab.  The  first  investigates  the  effects  of  increas¬ 
ing  A2e.  The  second  explores  the  possibility  of  improving 
the  model  fit  by  increasing  A23  and  Afe. 

Inspection  of  the  parameter  sets  shows  that  in  parameter  set 
1  (used  for  the  smaller  slab)  the  area  through  which  heat  is 
transferred  directly  from  the  slab  core  to  the  earth  (/t16) 
is  nearly  half  the  area  through  which  heat  is  transferred 
directly  from  the  slab  edge  to  the  earth  (A2g)  .  In  parame¬ 
ter  set  2,  however,  the  trend  is  reversed  so  that  Alb  is 
nearly  three  times  A2b.  While  it  is  anticipated  that  the 
effect  of  the  core  area  would  increase  as  the  slab  becomes 
larger,  this  dramatic  change  seems  excessive,  a  notion  which 
is  reinforced  by  the  apparent  underprediction  of  the  edge 
effect  from  parameter  set  1.  A  possible  source  of  error  is 
the  determination  of  the  distance  the  edge  effect  extends 
beyond  the  slab  edge  (d2) .  It  is  interesting  to  note  that 
this  distance  is  equal  to  the  characteristic  length  (-)  of 
the  smaller  slab.  It  is  postulated,  therefore,  that  d2 
should  be  defined  by  the  slab  characteristic  length  rather 
than  by  a  constant.  A  new  set  of  parameters  (see  Run  2B  in 
Table  3)  is  established  based  on  this  theory.  The  daily 
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average  flux  calculated  by  this  new  set  of  parameters  and 
the  difference  from  the  FDM  are  shown  in  Figure  15.  The 
curve  shows  a  better  fit  of  the  data  particularly  during  the 
winter  season  when  the  flux  is  highest.  Examination  of  the 
difference  between  the  FDM  and  GTF  Run  2B,  however,  shows 
that  the  model  still  demonstrates  a  systematic  seasonal 
error.  Comparing  Figure  14  and  Figure  15,  it  can  be  seen 
that  increasing  A2a  has  produced  a  linear  increase  in  the 
heat  flux  throughout  the  annual  cycle.  Nonetheless,  the 
model  is  reasonably  accurate  over  the  annual  cycle  giving  an 
error  of  3.6%  in  total  annual  energy  loss.  Table  4  gives  a 
numerical  comparison  of  GTF  models  2A  and  2B  with  the  base 
case  FDM. 

Table  4:  Results  of  Runs  2A  and  2B 
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ENERGY 

CONSUMP 

[%] 

TOTAL 
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ENERGY 
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FDM 

2.50 

44167.7 

RUN  2  A 

2.03 
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CO 
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Figure  15:  Flux  and  Difference  —  Run  2B 
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A /:i  and  A  fe  are  geometrically  similar  (see  eqns  51  and  52). 

The  correction  is  made  originally  on  A23  because  the  T2- 
A23-T3  system  is  more  stable  than  the  Te-Afe-Tf  system.  The 
correction  is  then  applied  to  both  A23  and  Afe.  The 
determination  of  the  corrected  A 23  is  based  on  the  assump¬ 
tion  of  a  functional  relationship  between  A23  and  the  slab 
characteristic  length.  Characteristic  length  is  discussed 
in  more  detail  in  Section  6.  v423  is  calculated  from  the 

equation 


A  ?3  2ur  e,ltllvL2d  (68) 

Assuming  that  A?/J  =  891m2  is  correct  for  the  small  slab 

A23  =  89l  =  ?.ur  oquivL?d  =  2nr  equiv(  1  1)  (69) 
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so  that 


891  (70) 

requw  2.U(ll)  J 

If  '  eq,nv  is  assumed  to  be  a  simple  function  of  characteristic 
length,  a  functional  relationship  can  be  postulated.  Sup¬ 
pose 

r~  equw=  K*  (characteristic.  Length)  (71) 

The  characteristic  length  (A/P)  of  the  small  slab  is  3.0m. 
So,  from  Equations  (70)  and  (71) , 

r  equiv  =  13=K*3.0  (72) 

so  that  K  is  approximately  equal  to  4.5.  Then,  using 
Equation  (71) ,  Equation  (68)  becomes 

3  2.1  =  l  ?  d(4  .l5)  (characteristic  Length).  (73) 

Similarly, 

A e/ =  2n  l.?e(  \.L3) (characteristic  Length).  (74) 

The  parameter  set  for  Run  2C  is  given  in  Table  3.  GTFs  and 
daily  average  heat  fluxes  are  recalculated  with  the  new  Az 3 
and  Aef.  The  numerical  comparison  of  Run  2C  to  the  FDM  and 
Runs  2 A  and  2B  is  tabulated  in  Table  5.  Run  2C  shows  a 
significantly  improved  fit  to  the  FDM  data  as  evidenced  by 
the  substantial  decrease  in  the  sum  of  squared  difference 
and  the  percentage  of  the  GTF  data  within  15%  of  the  FDM 
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data.  Unfortunately,  this  improvement  in  fit  was  accompa¬ 
nied  by  an  increase  in  the  total  error  over  the  annual 
cycle.  Examination  of  the  plot  of  the  difference  between 
the  results  of  RUN  2C  and  the  base  case  FDM  (Figure  16) 
shows  that  the  wintertime  droop  has  been  corrected,  but 
there  is  still  a  linear  underprediction  of  slab  loss. 


Figure  16:  Flux  and  Difference  —  Run  2C 
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It  was  demonstrated  by  Run  2B  that  such  an  error  might  be 
corrected  by  an  adjustment  of  A2e,  however,  an  adjustment  of 
the  magnitude  of  2B  is  unnecessary.  Once  again,  a  possible 
source  of  error  is  the  determination  of  the  distance  that 
the  edge  effect  extends  beyond  the  slab  edge  (d2) .  If  d2  is 
assumed  to  be  a  function  of  slab  characteristic  length  of 
the  form 


d?  ~  c  +  in (c lid r actor istic  length.) 


(75) 
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then  c  and  m  can  be  chosen  so  that  d2  varies  with  slab 

geometry.  Based  on  the  results  of  the  small  slab,  c  is 
probably  less  than  3.0  and  m  less  than  1.0.  If  c  is  assumed 
to  be  2 . 5m  then  m  can  be  calculated  from  the  model  of  the 
12m  x  12m  slab. 


which  gives 


3.0  =  2.5+  m(3.0) 


m  = 


3. 0-2. 5 

575 


(76) 


(77) 


which  is  approximately  equal  to  0.15.  When  the  equation 

d2  =  2.5  +  0. \l5(characteris’.lc  length)  (78) 

is  applied  the  new  values  of  d2  become: 
for  the  smaller  slab 


d?  =  2.9l5  (79) 

and  for  the  larger  slab 

d  2  =  4 .  1 9  (80) 

The  parameter  set  is  recalculated  and  given  in  Table  3,  Run 
ID.  When  this  change  is  made  in  the  calculations  for  the 
larger  slab  the  result  is  a  significantly  improved  fit  to 
the  FDM  base  case.  The  flux  and  difference  curves  are  shown 
in  Figure  17.  As  expected,  the  difference  curve  has  shifted 
linearly,  moving  nearer  to  the  zero  difference  line.  The 
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numerical  comparisons  are  given  in  Table  5.  The  quality  of 
the  fit  (as  measured  by  the  sum  of  squared  difference  and 
percentage  of  data  within  15%  of  the  FDM)  is  improved  as  is 
the  total  energy  consumption. 


Figure  17:  Flux  and  Difference  —  Run  2D 


Table  5:  Results  of  Runs  2A  -  2D 


MODEL 

MEAN 

FLUX 

[W/m'2] 

SUM  OF 

SQUARED 

DIFF 

%  OF 

DATA 

WITHIN 

15%  OF 

FDM  [%] 

TOTAL 

ANNUAL 

ENERGY 

CONSUMP 

%  ERROR 

IN 

TOTAL 

ENERGY 

CONSUMP 

[%] 

TOTAL 

ANNUAL 

DIFF  INS 

ENERGY 

CONSUMP 

FDM 

2.50 

44167.7 

RUN  2  A 

2.03 

106.9 

33 

35911.3 

-18.7 

-8256.4 

RUN  2B 

2.59 

54.0 

64 

45762.0 

+  3.6 

+1594.3 

RUN  2C 

19.9 

83 

40780.2 

-7.7 

-3387.5 

RUN  2D 

2.42 

8.1 

97 

42872.2 

-2.9 

-1295.5 

4.5.3  CORRECTED  MODEL  APPLIED  TO  12m  X  12m  SLAB 

Because  the  modifications  have  changed  the  functional  rela¬ 
tionships  which  effect  the  original  parameter  set,  it  is 
necessary  to  recalculate  the  parameter  set  for  the  smaller 
slab  and  recalculate  the  daily  average  fluxes  as  Run  ID.  A 
comparison  of  the  parameter  sets  of  Runs  1A-1D  is  given  in 
Table  6  and  the  results  in  Table  7.  The  quality  of  the  fit 
as  determined  by  the  sum  of  squared  difference  is  slightly 
degraded  by  the  changes  made  for  Run  ID.  Nonetheless,  the 
improvement  of  the  fit  of  the  model  to  the  data  for  the 
larger  slab  more  than  justifies  the  change. 
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Table  6:  Parameter  Sets  for  Runs  1A  -  ID 


RUN  1A 

RUN  IB 

RUN  1C 

RUN  ID 

AREA  [m2 J 

144 

144 

144 

144 

PERIMETER  [m] 

48 

48 

48 

48 

A/P  [m] 

3.0 

3 . 0 

3.0 

3.0 

157 

157 

297 

297 

^  12 

493 

493 

493 

493 

•*r 

CO 

r 

2S4 

284 

297 

23 

891 

891 

891 

933 

1  6  =  ^  1  rt 

105 

105 

105 

105 

2d  =  ^  2<> 

195 

195 

195 

192 

3  /  =  ^  3d 

867 

867 

867 

870 

l'l 

758 

758 

758 

758 

I7  2 

1416 

195 

195 

192 

l7  3 

6284 

867 

87 

87 

Table  7  Results  of  Runs  1A  -  ID 


MODEL 

MEAN 

FLUX 

[W/JIT2] 

SUM  OF 

SQUARED 

DIFF 

%  OF 

DATA 

WITHIN 

15%  OF 

FDM  [%] 

TOTAL 

ANNUAL 

ENERGY 

CONSUMP 

%  ERROR 

IN 

TOTAL 

ENERGY 

CONSUMP 

[%] 

TOTAL 

ANNUAL 

DIFF  IN 

ENERGY 

CONSUMP 

FDM 

6.13 

7716.9 

RUN  1A 

8.10 

2149.2 

37 

10194.4 

24.3 

+2477.5 

RUN  IB 

6.35 

214.2 

64 

7984.4 

MEM 

+267.5 

RUN  1C 

6.07 

69.9 

90 

7641.4 

-1.0 

-75.5 

RUN  ID 

6.19 

82.0 

89 

7786.6 

+0.9 

+69.1  | 
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5  FINAL  DEFINITION  AND  TESTING 


5.1  FINAL  DEFINITION 

The  model  used  for  Runs  ID  and  2D  gives  excellent  results 
for  both  the  12m  x  12m  square  slab  and  the  45m  x  45m  square 
slab.  The  parameter  sets  generated  for  both  these  runs  were 
constructed  from  the  same  series  of  equations  developed 
using  both  geometric  and  empirical  methods.  The  equations 
are  repeated  here  in  their  final  forms. 

D,  =  4.0m  -  0.5  rn  =  3.5  m.  (81) 

D?  =  15/7?  -  4/?7  =  llm  (82) 


D3  =  13.5m 

(83) 

(l  ,  =  1.0/71 

(84) 

d2  =  2.5  +  0. 1 5  (characteristic  length.) 

(85) 

DX*D2  „  ^ 

(86) 

h  = - -  =  7.25m. 

2 

ll2=  1  .0/77 

(87) 

li3  =  0. 1  m 

(88) 

1  lb  =  ~  !  3/  =  D  ) 

(89) 

II 

.V 

II 

II 

M 

(90) 

^  fe  t-  23  ^  3. 

(91) 
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(92) 


(93) 


^6,= 


2n(D,)(d,) 

.  r" 

In - - 


(94) 


Aef  =  2nA2e(4.5)  (characteristic  length ) 
23  i(D2)(d,) 


(95) 

(96) 


/1 23  =  2  u  Z.2rf(  4.5  )  (char  act  er  istic  lingth ) 

/^lb“/l]d  =  n(r  a_c^i)Z 

^  2e  =  ^  2d  =  Jt(r  a  +  ^2  )  _^16 

^3/  =  ^3d  =  ]l(ra  +  D3)  ~Aib~A2a 

^  ^  ,  a  , 

^2  “  ^Zl  h2 
'  3  ”  4  3  /  /*  3 


(97) 

(98) 

(99) 

(100) 
(101) 
(102) 
(103) 
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These  equations  are  used  to  construct  the  matrices 


A  = 


/ 


0 

A  12 
0 

0 


V 


A  l2  0  A  lb  0  A  ld 
0  A  23  0  0  A  2  d 

A  23  o  0  A3{  ■A  3d 

0  0  0  0  0 

0  A  3/  0  0  0 

A  2d  ^3d  0  0  0 

A2e  0  Abe  A  fe  0 


l 


/  = 


'■  12 
0 
A 


1  6 
0 


V1 


/ 


i  <( 


0 


L 12  0  Z.!b  0  /.  ld 

0  /,  23  0  0  L  2d 

L  23  ^  0  ^3/  ^  3d 

0  0  0  0  0 

0  A3/  0  0  0 

L?.d  1-3  d  0  0  0 

o  Abe  A  /e  0 


^  2e 
0 

^  be 

he 

0 

0 


l  = 


(104) 


(105) 


(106) 


The  geometry  matrices  A,  L  and  V,  along  with  the  soil 
property  matrices  k,  p,  and  cp  are  used  to  calculate  the 
matrices 


6' ,2  0 
0  6'  23 

£  23  0 
0  0 
0  c3/ 

^2d  ^  3d 

C2.  0 


k  lb 
0 
0 
0 
0 
0 


0  Cld  0  \ 

0  C2d  G  2e 

G  3  /  ^  3d  0 

0  0  cbe 

0  0  Gfe 

0  0  0 

0  0  / 


(107) 
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(108) 


by  the  equations 


(109) 


and 

C^ptcpiVu  <110> 

These,  in  turn,  are  used  to  generate  the  coefficient 
matrices 


(111) 
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(113) 


D  = 


0 

0 

G  be 


0 

~  G  3/  ~  G e/ 
0 

Gel 


0  C6e  \  (114) 

o  G„ 

~  G  ld  —  G Zd~  G 3d  0 

0  “  G  bt- G  t/  -  GZeJ 


which  are  used  with  Seem's  method  to  calculate  the  final 
multiple  input  GTFs  and  scalar  constants  which  for  the 
12m  x  12m  slab  are: 


2  0458  127  1  H9a*  000 
9  HI  I  1895231a  *002 
I  .7728083 1 44a ♦ 006 


-2. 184)0 /7414a*  00  7 
S.3595 1  80348a  *  008 
2.0349408  l  I0a*006 


5.3895 1 80348a*  005 
-8.9880223736a*  006 
1.3598085940a*  004 


2.05494081  I0a*006  \ 
1.3598085940a*  004 
-8. 5302667268a*  00  6/ 


(115) 


I  .  I  6985 l 3689 e  ♦  00/ 
2.ioor»076393c  ♦  001 
6 .317834 1 430c  •  00 
4.76281)  l 3469  c  ♦  00 


2 . 1 0050/6402 e  ♦  00  I 
6. 149 1664962ft +007 
-  5.41268197 1  I e ♦ 005 
5.505492RR56C  ♦  006 


-  6.3 1 7834  1 430  c  +  002 
-8.412681971  I  c  *  005 

2.4444345580c ♦ 007 

-  7.2500445257  e  ♦  003 


-  4.76285  I  3469c ♦  006  \ 

-  5.5054928856c+  006  \ 

-  7. 2500445257 e+ 003  I 
2.2994803858c ♦ 007  ) 


(116) 


f  1 .04 1 2722758ft +007 
I  .958  I  095233c  *  001 
9.6762874787c*  002 
\  4.2374081459c*  006 


l  .958 1 09S225c ♦ 00  I 
-5.76224  18860c +  007 
-  4.463  1067827c +005 
4.89687961  15c  +  006 


-  9.6762874788c ♦ 002 

-  4.463  1 067827  c  +  005 

-  2.204993 1343c  +  007 

-  1. 4768777916c*  004 


4.2374081459c* 006  \ 
4.89687961  lSc  +  006  | 

-  1.4768777916c* 004  I 
-2.0523358176c*  007/ 


(117) 


/  3.0658049564c  +  006 
I  .6560367 454c  + 000 
6.77 15988886c*  002 
V-  I  247  I  449406c*  006 


1 .6560367478c*  000 
1 .797  1 572628ft ♦ 007 
4.52  I  I  655346a  *  005 
-  I  4460026224e*006 


6.77  15988886a  *002 
4.521  1655346  8*005 
6.59295161018*006 
8.5289  I  48747  a  *003 


-  1.247  I  449406a*  006  \ 
-!. 4460026224a* 006  | 
8.5289  148747a*  003  I 
6.0581677786e*006  ) 


(118) 


c,=  -2.686/4  (119) 

p2  =  2.39032 
e-,  =  -0. 703485. 
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5.2  VALIDATION 


The  final  set  of  GTF  coefficients  and  scalar  constants  cal¬ 
culated  using  the  above  equations  are  used  to  test  the  model 
for  a  variety  of  conditions  including  diverse  climates,  slab 
size  and  shape,  and  sensitivity  to  input  data. 

5.2.1  8IZE 


The  effect  of  slab  size  on  the  accuracy  of  the  model  has 
been  described  in  Section  4.5.  In  summary,  the  model  is 
quite  accurate  for  relatively  small  (144  m2)  to  relatively 
large  (2025  m2)  square  slabs  giving  an  error  in  total 
annual  energy  consumption  of  less  than  3%  in  both  cases. 
The  model  is  slightly  more  accurate  overall  for  the  larger 
slab  based  on  the  percentage  of  the  data  within  15%  of  the 
FDM:  97%  for  the  larger  slab  vs.  89%  for  the  smaller  slab. 


5.2.2  CLIMATE 

The  final  GTF  coefficients  and  scalar  constants  were  used 
with  environmental  data  for  Medford  OR,  Philadelphia  PA,  and 
Phoenix  AZ.  Plots  of  the  flux  and  difference  for  all  four 
locations  are  given  in  Figures  18,19,20  and  21.  Table  8 
gives  numerical  data  regarding  the  accuracy  of  the  models. 
For  Minneapolis,  Medford,  and  Philadelphia,  the  difference 
between  the  GTF  model  and  the  FDM  is  very  nearly  zero.  In 
all  cases  the  difference  is  less  than  1  UWm2  except  for  a 
few  days  at  the  beginning  of  the  annual  cvcle.  In  Phoenix 
where  the  annual  mean  flux  is  approximately  1.5  W/mz,  an 
error  of  less  than  1  W  /  m2  can  create  a  significant  error 
when  the  actual  value  of  the  error  is  quite  small. 
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Figure  18:  Flux  and  Difference  —  Minneapolis  MN 


Figure  19:  Flux  and  Difference  —  Medford  OR 
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Figure  20:  Flux  and  Difference 


Philadelphia  PA 


Figure  21:  Flux  and  Difference  —  Phoenix  AZ 


Table  8:  Results  of  GTF  Models  for  Various  Locations 


MODEL 

LOCATION 

MEAN 

FLUX 

[W/m*  2 

%  OF 

DATA 

WITHIN 

15%  OF 

FDM 

[%] 

TOTAL 

ANNUAL 

ENERGY 

CONSUMP 

[kWhr ] 

%  ERROR 

IN 

TOTAL 

ENERGY 

CONSUMP 

[%] 

TOTAL 

ANNUAL 

DIFF  IN 

ENERGY 

CONSUMP 

[ kWhr ] 

FDM 

Minneapolis 

6.10 

7716.9 

GTF 

Minneapolis 

6.19 

89 

7786.6 

+  0.9 

+  69.7 

FDM 

Medford 

3.87 

4867 . 1 

GTF 

Medford 

3 .86 

78 

4856.6 

-0.2 

-10.5 

FDM 

Philadelphia 

3 . 89 

4893 . 7 

GTF 

Philadelphia 

3.87 

78 

4863.0 

-0.6 

-30.7 

FDM 

Phoenix 

-1.66 

-2082.9 

GTF 

Phoenix 

mm 

72 

-1832.7 

-12.0 

+250.2 

5.2.3  SHAPE 

This  model  was  developed  assuming  a  square  slab  and  uses  the 
circular  isotherms  which  evolve  as  the  result  of  that  geome¬ 
try.  Although  it  was  not  expected  that  this  model  would 
adequately  model  non-square  slabs,  the  extent  of  the 
inaccuracy  was  unknown.  Therefore,  parameter  sets  were  con¬ 
structed  based  on  the  slab  perimeters  and  areas  and  using 
the  above  equations.  These  parameter  sets  were  used  to 
calculate  GTF  coefficient  matrices  and  scalar  constants,  and 
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from  them,  daily  average  heat  fluxes.  The  results  are  shown 
in  Figures  22  and  23.  Numerical  comparison  of  these  results 
to  the  FDM  results  for  the  non-square  slabs  is  given  in 
Table  9. 

As  was  anticipated,  the  model  does  not  give  good  results  for 
non-square  slabs.  The  form  of  the  errors  indicates  an  inac¬ 
curacy  in  the  calculation  of  the  edge  effect.  In  the  summer 
when  the  area  effect  dominates,  the  difference  between  the 
FDM  and  the  GTF  model  is  nearly  zero.  However,  as  the 
ground  surface  temperature  drops  and  the  edge  effect  becomes 
more  important,  the  difference  between  the  FDM  and  the  GTF 
model  shows  an  increasing  underprediction  of  the  slab  heat 
loss . 


Figure  22:  Flux  and  Difference  —  6  x  24  Rectangle 
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Figure  23:  Flux  and  Difference  —  18  x  112  Rectangle 
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Table  9:  Results  of  GTF  Model  for  Non-square  Slabs 


MODEL 

SLAB  SIZE 

[itT2] 

MEAN 

FLUX 

[W/m‘2 

FDM 

6  x  24 

7 . 30 

GTF 

6  x  24 

5.81 

FDM 

18  x  112 

3.19 

GTF 

18  X  112 

2.18 

FLUX  DATA  ANNUAL  IN  ANNUAL 

[W/m‘2  WITHIN  ENERGY  TOTAL  DIFF  IN 

15%  OF  CONSUMP  ENERGY  ENERGY 

FDM  [%]  [kWhr]  CONSUMP  CONSUMP 

[%]  [ kWhr ] 


9180.2 


7305.2  -25.7  -1875.0 


56405.2 


38530.3 


-17874.9 
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5.2.4  SENSITIVITY  TO  INPUTS 


It  is  important  to  understand  the  effect  of  the  accuracy  of 
the  input  data  on  the  results  of  the  model  particularly  if 
the  required  data  are  not  available  and  approximations  must 
be  made.  The  most  probable  approximations  are: 

(1)  substituting  daily  average  outdoor  dry  bulb  tempera¬ 
ture  for  daily  average  ground  surface  temperature  (T ,) 

(2)  substituting  annual  average  outdoor  dry  bulb  temper¬ 
ature  for  annual  average  ground  surface  temperature  (T d) 

(3)  substituting  constant  indoor  air  temperature  for 
daily  average  floor  surface  temperature  (T b)  and  daily 
average  floor  edge  temperature  (7'p)  . 

These  approximations  will  be  tested  for  the  12m  x  12m  slab 
with  the  four  climatic  conditions  and  the  45m  x  45m  slab 
with  Minneapolis  MN  climate. 

The  following  runs  used  the  final  GTF  coefficients  and  sca¬ 
lar  constants.  When  the  approximation  for  one  input  data 
set  was  used,  the  remaining  inputs  were  held  identical  to 
those  in  Runs  ID  and  2D. 

GROUND  SURFACE  TEMPERATURE 

Table  10  gives  the  numerical  comparison  of  the  data  result¬ 
ing  from  substituting  the  daily  average  outdoor  air  tempera¬ 
ture  for  daily  average  ground  surface  temperature  as  the 
input  at  7/.  The  error  in  total  energy  consumption  over  the 
entire  cycle  ranges  from  17.1  %  to  219.5  %.  Inspection  of 
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the  graphical  representation  of  the  data  (Figures  24,  25, 

26,  27  and  28)  reveals  a  common  pattern  in  the  error.  In 
all  cases,  the  largest  factor  in  the  error  is  a  positive 
linear  offset  which  is  greatest  in  Phoenix  where  the  temper¬ 
ature  difference  between  the  air  temperature  and  ground  tem¬ 
perature  is  highest.  The  larger  slab,  where  the  edge  effect 
is  less  substantial,  shows  a  much  smaller  effect  of  changing 
far  field  temperature. 

DEEP  GROUND  TEMPERATURE 


Table  11  and  figures  29,  30,  31,  32  and  33  show  the  results 
of  using  the  annual  average  outdoor  air  temperature  as  the 
deep  ground  temperature,  T d.  Because  the  annual  mean  out¬ 
door  air  temperature  is  less  than  the  annual  mean  ground 
surface  temperature,  the  calculated  heat  loss  from  the  slab 
is  correspondingly  higher.  Again,  the  largest  component  of 
the  error  is  the  linear  offset.  The  error  is  small  in 
moderate  and  cold  climates  —  less  than  5%.  In  Phoenix 
there  is  only  a  slight  difference  between  the  annual  mean 
outdoor  air  temperature  and  the  indoor  air  temperature 
( 0 . 1C)  so  that  the  anrual  mean  heat  flux  through  the  slab  is 
very  small  compared  to  other  locations  where  the  temperature 
difference  is  an  order  of  magnitude  larger.  At  the  same 
time,  the  difference  between  the  annual  mean  outdoor  air 
temperature  and  the  annual  mean  ground  surface  is  much 
greater  in  Phoenix  due  for  the  most  part  to  the  increased 
solar  gain.  The  consequence  of  these  two  conditions  is  a 
much  greater  effect  from  this  substitution  than  is  seen  in 
the  other  climates.  In  general,  however,  it  gives  good 
results  and  even  gives  a  slightly  better  fit  in  some  cases 
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as  evidenced  by  the  increase  in  the  percentage  of  the  data 
within  15%  of  the  FDM  data  for  small  slab  in  Medford  and  the 
larger  slab  in  Minneapolis. 

SLAB  TEMPERATURE 

In  this  case  a  constant  value  is  substituted  for  the  input 
data  set.  This  is  a  convenient  substitution  and  practical 
for  the  many  cases  where  the  slab  temperature  is,  in  fact, 
nearly  constant.  Table  12  and  Figures  34,  35,  36  and  37 
give  the  results  of  this  substitution.  In  assuming  a  con¬ 
stant  temperature  approximately  10%  higher  than  the  actual 
floor  surface  temperature  a  error  of  roughly  10%  is 
introduced.  This  error  is  primarily  a  linear  shift,  which 
appears  typical  of  input  data  set  changes.  It  seems  to  be 
due  for  tne  most  part  to  the  difference  between  the  mean 
value  of  the  original  data  set  and  the  mean  of  the  substi¬ 
tuted  data.  As  in  all  the  other  cases  of  input  substitu¬ 
tion,  the  effect  is  substantially  smaller  for  the  larger 
s’  b. 

GENERAL  COMMENTS 


For  the  most  part,  the  changes  described  in  this  section 
cause  linear  shifts  of  the  flux  curve.  This  linear  shift 
appears  to  be  related  principally  to  the  difference  between 
the  mean  of  the  original  data  set  and  the  mean  of  the  sub¬ 
stituted  data.  The  slight  changes  in  the  shape  of  the  input 
data  curves  do  not  have  a  great  effect  on  the  flux.  It  is 
probable,  based  on  the  linearity  of  the  change,  that  alter¬ 
ing  more  than  one  input  would  result  in  a  linear  shift 
related  to  the  added  effects  of  the  individual  changes. 
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Table  10:  Results  of  Substituting  Daily  Average  Outdoor  Air 
Temperature  for  Daily  Average  Ground  Surface  Temperature 


MODEL 

LOCAT ' N 

EDGE 

SIZE 

[m] 

MEAN 

FLUX 

[W/m'2] 

%  OF 

DATA 

WITHIN 

15%  OF 

FDM 

[%] 

TOTAL 

ANNUAL 

ENERGY 

CONSUMP 

[ kWhr ] 

%  ERROR 

IN 

TOTAL 

ENERGY 

CONSUMP 

m 

TOTAL 

ANNUAL 

DIFF  IN 

ENERGY 

CONSUMP 

[ kWhr] 

EDM 

Minn 

12 

D 

7716.9 

GTE 

Minn 

12 

8.11 

mm 

10197.4 

+  24.3 

+2480.5 

FDM 

Medford 

12 

3.87 

4867 . 1 

GTF 

Medford 

12 

6.3  0 

1  7 

7923 . 0 

+  62.8 

+  3055.9 

FDM 

Phi  la 

12 

3 . 89 

4893.7 

GTF 

Phi  la 

12 

5.65 

12 

7101.6 

+  45.1 

+2207.9 

FDM 

Phoenix 

12 

-1.66 

-2082.9 

GTF 

Phoenix 

□ 

1.98 

2 

2488.8 

-219.5 

+4571.7 

FDM 

Minn 

m 

2.50 

44167.7 

GTF 

Minn 

m 

2.92 

47 

51711.8 

+  17.1 

+7544 . 1 
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Figure  24:  Flux  and  Difference  —  —  Minneapolis  MN 
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Figure  25:  Flux  and  Difference  —  I  ,  -  / —  Medford  OR 
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Figure  26:  Flux  and  Difference  —  / -  7 —  Philadelphia 

PA 
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Figure  27:  Flux  and  Difference  —  7  /  /,,,,  —  Phoenix  AZ 
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Figure  28:  Flux  and  Dif  ference  —  I  /  -  1  45  x 

Minneapolis  MN 
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Table  11:  Results  of  Substituting  Annual  Average  Air  Tem¬ 
perature  for  Annual  Average  Ground  Surface  Temperature 


modi:  i. 

LOCAT'N 

EDGE 

SIZE 

[m] 

MEAN 

FLUX 

[ W/m ‘ 2 ] 

%  OF 

DATA 

WITHIN 

1  Si  OF 

FDM 

m 

TOTAL 

ANNUAL 

ENERGY 

CONSUME 

[kWhr] 

%  ERROR 

IN 

TOTAL 

ENERGY 

CONSUMP 

[%) 

TOTAL 

ANNUAL 

DIFF  IN 

ENERGY 

CONSUMP 

f  kWhr ] 

FDM 

GTF 

FDM 

Minn 

Minn 

12 

6 .  1  3 

7716.9 

1  2 

<>.  12 

87 

7  9 SO  . 

12.9 

t  23  3  .  S 

Med  f  ord 

1  2 

3  .  87 

4867 . 1 

GTF 

Med  f  ord 

12 

c 

79 

SOS  3 . 7 

1  4  .  S 

3216.6 

FI)M 

Phi  1  a 

12 

3 . 80 

4893.7 

GTF 

I’h  i  1  a 

12 

3.99 

7  8 

S02  0 . 4 

warn 

4126.7 

FI)M 

Phoen i x 

12 

-1  .  66 

-2082 . 9 

GTF 

FDM 

Phoen i x 

12 

-1 . 03 

12 

-1299.7 

-37 . 6 

4  783. 2 

Minn 

1 

2  .  SO 

44167.7 

GTF 

M  i  nn 

4S 

2 . 5 1 

100 

4  4  3  9  S  .  6 

+  0. 5 

+  227.9 
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Figure  31:  Flux  and  Difference  —  =  Annual  Mean  T 
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Figure  32:  Flux  and  Difference  —  =  Annual  Mean  / 
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Figure  33:  Flux  and  Difference  —  T d  =  Annual  Mean  7 

45  x  45  —  Minneapolis  MN 
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Table  12:  Results  of  Substituting  Constant  Indoor  Air  Tem¬ 
perature  for  Daily  Average  Slab  Surface  Temperature 


MODEL 


FDM 

GTF 

FDM 

GTF 

FDM 

GTF 

FDM 

GTF 

FDM 

GTF 


Figure  34:  Flux  and  Difference  —  Tb 
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Figure  35:  Flux  and  Difference  —  Tb  =  Tia  —  Medford  OR 
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Figure  36:  Flux  and  Difference  —  Tb  =  Tia  —  Philadelphia  PA 


Figure  37:  Flux  and  Difference  —  T b  =  T  ,a  —  Phoenix  AZ 
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Figure  38:  Flux  and  Difference  — 

Minneapolis  MN 


6  NETWORK  PARAMETERS  BASED  ON  CHARACTERISTIC  LENGTH 

Bahnfleth  [3]  reports  that  heat  flux  through  slabs  of  sev¬ 
eral  different  rectangular  geometries  can  be  calculated 
based  on  the  slab  characteristic  length  (A/P) .  This  leads 
to  the  speculation  that  it  may  be  possible  to  define  the 
network  parameters  as  functions  of  the  soil  geometry  and 
slab  characteristic  length  (A/P) ,  thereby  allowing  the  use 
of  the  model  with  non-square  slabs. 

As  a  crude  test  of  this  proposition,  empirical  models  were 
developed  for  slabs  of  four  different  configurations:  12m  x 
12m,  45m  x  45m,  6m  x  24m,  and  18m  x  112m.  Little  attempt 
was  made  at  this  point  to  attach  geometric  significance  to 
the  network  parameters,  but  rather,  each  parameter  set  was 
adjusted  based  primarily  on  the  quality  of  the  resulting  fit 
to  each  individual  set  of  base  case  data.  A  parameter  set 
was  considered  acceptable  when  more  than  80%  of  the  result¬ 
ing  data  were  within  15%  of  the  corresponding  FDM  data  and 
the  error  in  total  annual  heat  flux  was  less  than  5%  with 
the  Minneapolis  MN  weather  data.  It  should  not  be  assumed 
that  these  parameter  sets  are  in  any  way  optimal.  Once  a 
set  of  parameters  for  each  configuration  was  developed,  the 
network  parameters  were  compared  in  an  effort  to  identify 
patterns  among  the  four  cases  (see  Table  13) . 

Several  relationships  became  evident.  /13/.  A?3,  and  \  2 

are  of  identical  or  similar  value  when  the  area  of  the  slabs 
are  (nearly)  identical.  This  is  a  strong  indication  that 
those  parameters  are  functions  of  the  slab  area.  Corre¬ 
spondingly,  .l(,„  and  \2d  appear  to  be  functions  of  the  slab 
perimeter.  The  remaining  parameters  are  assumed  to  be 
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functions  of  characteristic  length,  A/P.  Pairing  the  two 
square  slabs  and  the  two  non-square  slabs,  a  line  is  fit  to 
each  set  and  the  coefficients  of  the  resulting  equations 
compared.  For  equations  of  the  form 

y  =  ci  +  bx  (114) 

the  coefficients  and  variables  are  shown  in  Table  14. 

It  is  clear  that  a  single  set  of  linear  equations  for  the 
network  parameters  in  terms  of  the  slab  area,  perimeter  or 
characteristic  length  can  be  written  and  should  give  accept¬ 
able  results  for  all  four  cases.  A  suggested  set  of  equa¬ 
tions,  based  on  these  data  is  shown  in  Table  15.  Network 
parameters  are  calculated  using  these  data  and  are  shown  in 
Table  16. 

Using  the  original  L  matrix  and  the  A  and  V  matrices  gener¬ 
ated  using  Table  16,  new  GTFs  and  scalar  constants  are 
calculated  and  QCALC  is  used  to  calculate  the  daily  average 
heat  flux  through  the  slab  using  Minneapolis  MN  climatologi¬ 
cal  data.  The  results  are  shown  in  Figures  38,  39,  40,  and 
41.  Numerical  comparisons  are  given  in  Table  17.  In  all 
cased  more  than  80%  of  the  data  are  within  15%  of  the  FDM 
data,  and  the  error  in  total  energy  consumption  is  less  than 
10%.  It  is  evident  that  it  is  possible  to  develop  a  set  of 
equations  for  calculating  network  parameters  as  functions  of 
slab  area,  perimeter  and  characteristic  length  which  give 
good  results  for  a  variety  of  rectangular  geometries.  The 
equations  developed  for  this  example  should  not  be  consid¬ 
ered  universal;  a  more  rigorous  method  of  parameter  estima- 
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tion  should  be  used  to  develop  a  truly  generic  parameter 
set.  Nonetheless,  this  example  indicates  that  such  a 
procedure  should  yield  good  results. 

Table  13  Parameter  Sets  for  Empirical  Models 


AREA  [  m2] 


PERIMETER  [m] 


A/P  [m] 


12x12 


144 


48 


3.00 


87 


45x45 


2025 


180 


11.25 


95 


6x24 


144 


60 


2. 


82 


18x112 


2016 


260 


7.75 


97 


-0.6  -3.8 


157 

608 

121 

415 

493 

1934 

385 

1322 

200 

650 

255 

1000 

800 

3000 

792 

3000 

100 

1881 

90 

17*9 

400 

1500 

600 

2600 

800 

2000 

878 

2000 

80 

1000 

80 

1000 

320 

1850 

320 

2000 

100 

100 

100 

100 

80 


Table  14:  Equations  of  Lines  Fit  to  Empirical  Data 
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Table  15:  Equations  of  Common  Lines  Fit  to  Empirical  Data 


PARAMETER 

a 

b 

VARIABLE 

*/. 

-9.0 

55 

A/P 

AI2 

-33 

175 

A/P 

34 

3 . 6 

P 

A  23 

627 

1.2 

A 

A  |  6  =  ^  1  d 

-36 

.9 

A/P 

2d  ~  2 o 

0 

9 

P 

A  3/=  A  2d 

750 

.6 

A 

V  i 

9.4 

.5 

A 

^2 

203 

.9 

A 

^3 

100 

0 

CONSTANT 

82 


Table  16:  Parameter  Sets  Calculated  Using  Empirical 

Equations 


BHHHI 

12X12 

45x45 

6X24 

18x112 

AREA  [m2] 

144 

2025 

144 

2016 

PERIMETER 

[m] 

48 

180 

60 

260 

A/P  [m] 

3.00 

11.25 

2.4 

7.75 

A*. 

156 

610 

123 

417 

^  1  2 

492 

1936 

387 

1324 

'  ^  6e 

207 

682 

250 

970 

^  23 

800 

3057 

800 

3046 

mam 

94 

1787 

94 

1778 

2d  ^  2c 

432 

1620 

540 

2340 

'^3/  =  ^3d 

836 

1965 

836 

1960 

l/, 

81 

1022 

81 

1017 

1/2 

327 

2020 

327 

2011 

I  "' 3 

100 

100 

100 

100 
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Figure  39:  Flux  and  Difference  —  Run  3A 
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Figure  40:  Flux  and  Difference  —  Run  3B 
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Figure  41:  Flux  and  Difference  —  Run  3C 


Table  17:  Results  of  GTF  Model  Using  Parameter  Sets  Based 

on  Empirical  Equations 


"  - - - - - - =i 

MODEL 

MEAN 

FLUX 

[W/m'2] 

%  OF 

DATA 

WITHIN 

15%  OF 

FDM 

[%] 

. .  .  '  ■ 

TOTAL 

ANNUAL 

ENERGY 

CONSUMP 

[kWhr ] 

%  ERROR 

IN 

TOTAL 

ENERGY 

CONSUMP 

[%] 

TOTAL 

ANNUAL 

DIFF  IN 

ENERGY 

CONSUMP 

[ kWhr ] 

FDM- 12x12 

6.10 

7716.9 

GTF-12X12 

6.57 

80 

8263 . 1 

+  6.6 

+546.2 

FDM-6X24 

7.30 

9180.2 

GTF-6X24 

7.31 

81 

9190.5 

+  0.1 

+  10.3 

FDM-45X45 

2.50 

44167.7 

GTF-45X45 

2.47 

93 

43763.6 

-0.9 

-404 . 1 

FDM-18X112 

3.19 

56405.2 

GTF-18X112 

2.93 

86 

51765.1 

-8.2 

-4640.1 
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7  UTILIZATION  OF  THE  GTF  MODEL  FOR  ENERGY  ANALYSIS 

Because  of  its  conceptual  similarity  to  existing  energy 
analysis  programs  using  transfer  function  models  of  building 
components,  this  model  is  particularly  suitable  for  incorpo¬ 
ration  into  these  programs.  When  used  with  these  types  of 
hourly  energy  analysis  programs,  the  ground  network  would  be 
seen  as  another  zone  connected  to  the  conditioned  space  by 
an  "interzone  partition"  which  would  include  the  slab  itself 
and  the  top  0.5m  of  soil.  The  surface  inner  and  outer  tem¬ 
peratures  of  the  "partition"  would  be  the  calculated  hourly 
slab  surface  temperature  and  the  GTF  model  slab  center  and 
edge  temperatures  (which  are  equal  in  this  case) ,  respec¬ 
tively.  An  algorithm  for  the  calculation  of  ground  surface 
temperatures  (such  as  TEARTH,  developed  by  Bahnfleth)  must 
be  included  in  the  processing  of  weather  data  in  order  to 
provide  correct  input  values  for  the  far-field  and  deep 
ground  temperatures. 

This  model  could  also  serve  as  part  of  a  stand-alone  slab 
heat  loss  program  for  situations  where  the  daily  average 
slab  surface  temperatures  are  Known  or  can  be  reasonably 
approximated.  For  example,  Section  5.2.4  showed  that  sub¬ 
stituting  a  constant  indoor  air  temperature  in  place  of  the 
slab  surface  temperature  gave  acceptable  results  for  all 
locations  and  both  slab  sizes.  The  program  would  also 
require  an  algorithm  for  the  calculation  of  undisturbed 
ground  surface  temperatures. 
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8  CONCLUSIONS 


A  simple  multiple-input  transfer  function  model  of  the  heat 
transfer  in  the  ground  beneath  a  square  slab  was  presented. 
It  was  tested  and  modified  to  model  both  relatively  small 
slabs  where  edge  effects  are  strong  and  larger  slabs  whose 
heat  flux  is  more  strongly  effected  by  the  flux  through  the 
core.  Tested  over  a  broad  range  of  climatic  conditions,  the 
model  calculates  slab  heat  flux  within  1  W'  /m2  at  all  times 
and  for  all  locations.  This  translates  to  an  error  of  less 
than  1%  (as  compared  to  the  detailed  finite  difference 
model)  for  moderate  and  cold  climates  and  12%  for  Phoenix 
where  the  total  flux  is  very  low.  The  accuracy  of  the  model 
is  dependent  upon  the  accuracy  of  the  input  data,  however, 
some  reasonable  approximations  to  the  necessary  input  data 
can  give  acceptable  results. 

The  full  capability  of  the  model  was  not  tested  in  this 
study.  Further  work  to  develop  a  definition  of  the  network 
parameters  based  on  characteristic  length  could  expand  the 
use  of  the  model  to  non-square  and  possibly  even  non- 
rectangular  surfaces.  Testing,  and  if  necessary,  modifica¬ 
tion  of  the  parameter  equations  to  support  differential  slab 
core  and  slab  edge  temperatures  would  allow  the  model  to  be 
used  more  effectively  for  insulation  studies  where  the 
placement  of  the  insulation  is  contingent  on  the  slab  to 
environment  temperature  difference. 
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APPENDIX  A  SEEM'S  METHOD 


Seem’s  Method  for  Calculating  Transfer  Functions 
for  Multidimensional  Heat  Transfer 


The  system  is  described  by  a  series  of  equations  in  state 
space  formulation: 


a  A 

JT 


=  ,1  A  +  3U 


(A.l) 


Q  =  c  A  +  ©1/ 


(A. 2) 


where 


V  =  vector  of  state  variables 


£/ =  vector  of  known  input  variables 


Q=  vector  of  output  variables, 


The  solution  of  this  system  of  equations  is: 


,Y,.6»<?*6.Y,  +  J  oA""~"BV(x)dx 


•(-6 


-<  (  (  •  i  -  l  ) 


The  input  function  of  the  form: 


(A. 3) 


90 


uw-u, 


x  - 1 


(A. 4) 


6 


(-6 


*/.) 


is  substituted  into  the  solution  equation  and  the  integrals 
evaluated  leaving  the  solution  equation 


A 


4>*i  +  (r,-r2)i/,-*T2£/l*5 


(A. 5) 


where 


<t>  =  e*6  =  the  exponential  matrix 
T,  =  A~l  (e*b  -  l)S 


The  forward  shift  operator 

Fvt  =  iVs  (A. 6) 

which  relates  u,+6  to  its  previous  values  vt,  is  used  in  the 

solution  equation  in  order  to  \.sr  the  equation  entirely  in 
time  t: 


(f7-4>),Y,«(fT2  +  ri  -r2)U,  (A. 7) 

or 

.vl-(r/-4>)*,(rr2  +  r,-r2)i/l. 

This  equation  gives  the  state  variables  in  terms  of  the 
inputs.  It  is  substituted  back  into  the  state  output 
equation  to  give  the  outputs  in  terms  of  the  inputs: 
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Qt  =  (C(Fi-<t>yl(r\ ■,-rl-r2)  +  2>)u( 


(A. 8) 


where 


(FI 


R0Fn-'  +  R)F""  +  ...+  R„-2F+Rn.i 


■n-2 


Fn+e,F 


n- 1 


+  g„.  \  F  +  Q, 


(A. 9) 


Substituting,  combining  common  terms  of  the  forward  shift 
operator,  F,  and  shifting  the  equation  n  timesteps  back 
leaves  the  transfer  function  equation 


Q, 


(C  R0\~  ?_  +  <D)U  t 


(A. 10) 


+  (C(R0(i \  -  r2)+  /?,  r,) +  <vPU''/-6 

-(C(fl,(rt  -r2)  +  /??r2)  +  G2£> . .. 

+  (£(^n-?0  1  “  f  ?.)+  R  n  -  1  I  2  )  +  G  |  "D)  U  , .  (n  .  i  )b 
-(CRn_i(\\-V,)-onV)Ul.n,. 

More  concisely. 


Q,= 


£(s;fv, 


(A. 11) 


where 

s0  =  c/?0r '*  +  *> 

s;=-c(«;.,(r,-r7)-«/r?)-e/p  /on  </<n-  i 
5n  =  c/?n.,(r,  -r2)  +  Gn>r» 


Seem's  computationally  efficient  procedure  for  calculating 
numerically  significant  coefficients  uses  Leverrier's  algo¬ 
rithm  with  the  analytical  solution  as  follows: 
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l.  Calculate  the  exponential  matrix  by  a  truncated  power 
series  expansion  with  scaling  and  squaring 

1.1  Calculate  the  matrix  row  norm  of  matrix  .36 


|.36|«,=  max  £  |tti/6| 

1  £  t  S  n  /  “  1 


(A. 12) 


1.2  Find  the  smallest  integer  k  such  that 

2*>|.36|„  (A. 13) 

1.3  Divide  matrix  .36  by  2*. 

1.4  Calculate  the  number  of  terms  to  keep 

/.  =  min  {(3|.36|  +  6)v  100}  (A. 14) 

Ab 

1.5  Calculate  e2* 


1 

^  f  *  6  \  3 

MM 

.3  6  1 

_  4-  . 

4  4 

2* 

2! 

3! 

L\ 

(A. 15) 


1.6  Calculate  o'*6 

j,  t  (A'!6) 

2.  Calculate  T,  and  V ? 
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r,  =  .v'cg*6-/)® 


(A. 17) 


-i  >'i 

r2  =  .q  \-~-3) 
o 

3.  Calculate  the  S0  matrix 

S0  =  CR0T?  +  'D  (A.  18) 

4.  Calculate  the  S  matrices  for  1  < j  <  n  -  I  where  n  is  the 
number  of  state  variables 

4.1  Set  starting  R  matrix  to  the  identity  matrix 

R  =/  (A.  19) 

4.2  Calculate  the  scalar  constant,  e,  and  the  next  R 
matrix,  Rriou, 

Rou,  =  Rnma  (A*2°) 

Trace(4>KaW) 

G/  = - - - 

«„TO  =  <l,«„w  +  o// 

4.3  Calculate  the  next  S  matrix 

■S ,-C{  Rout ( I' ,  -  I'z )  +  r2 )  +  e / £>  (A.  21) 

4.4  Iterate  4.2  -  4.3  resetting  j  =  j+1  until  j  = 
n-1 
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5.  Calculate  scalar  constant  and  S  matrix  for  j 


T  r  a  co(<t> 

°"  = - - - 


n 

(A. 22) 
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APPENDIX  B  TRUE  BASIC  PROGRAM  GTF 


TRU BASIC  PROGRAM  GTF 

FOR  CALCULATING  MULTIPLE  INPUT  GROUND  TRANSFER  FUNCTION 
COEFFICIENTS  AND  SCALAR  CONSTANTS 

RECORD, RECSIZE  64 

INPUT  PROMPT  "RUN  ID?  runid$ 

OPEN  #3:  NAME  "D: \TRUBASIC\INDATA\"&runid$&"GTFs" , CREATE 
NEWOLD, ORGANIZATION 
SET  #3:  POINTER  BEGIN 
WRITE  #3:  DATE$ , TIME$ 

WRITE  #3:  runid$ 

DIM  CONDUCTIVITY (7, 7) , DENSITY (3) , SPECIFIC_HT ( 3 ) 

DIM  AREA (7,7), LENGTH (7,7) 

DIM  P ( 4 ) , H ( 7 ) , AOV ( 7 ) , VOLUME ( 3 ) 

DIM  DISTANCE (4) 

DIM  G ( 7 , 7 ) , CINV ( 3 ) ,D(4,4),C(4,3),A(3,3),B(3,4) 

DIM  I DEN (3,3) ,PHIADJ(3,3) , PHI (3, 3) ,PHINEW(3,3) , FACTOR (3, 3) 
DIM  ADEL (7,7), AADJ (7,7) 

DIM  GAMMA 1 (3,3), GAMMA 2 (3,3), AINV (3,3), TEMPI (4,4), TEMP2 (4,4) 
DIM  ROLD(10, 10) , RNEW ( 10 , 10) , SOLD (10, 10) ,SNEW(10, 10) ,SFINAL(1 
0,10) 

DIM  TABLE (4, 4) 

CALL  PROPERTIES 

CALL  MODEL 


REM  Generic  Equation  Construction 
REM  set  up  equation  matrices 
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REM  set  up  general  configuration 

REM  for  model  with  3  state  variables  and  4  inputs 
REM 

LET  num_state_var=3 
LET  num_inputs=4 


REM 

REM  Generate  thermal  resistance  matrix  G  with  equation 
REM  G(i,j)  =  k(i, j)*A(i, j)/L(i, j) 

REM 

MAT  G=ZER 

FOR  i=l  to  num_state_var  +  num_inputs 
FOR  j=l  to  num_state_var  +  num_inputs 
IF  LENGTH ( i , j ) =0  THEN 
LET  G(i,j)=0 
ELSE 

LET  G ( i , j ) =CONDUCTIVITY (  i ,  j  ) *AREA (  i  ,  j  ) /LENGTH ( i , j ) 
END  IF 
NEXT  j 
NEXT  i 

REM 

REM  Generate  the  matrix  of  the  inverses  of  the  thermal 
capacitances 

REM  CINV(i)  =  (DENSITY (i) *SPECIFIC_HT (i) *VOLUME(i) ) A (-1) 
REM 

FOR  i=l  to  num_state_var 

LET  CINV(i)= (DENSITY (i) *SPECIFIC_HT ( i) *VOLUME(i) ) * (-1) 

NEXT  i 
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REM 

REM  Generate  the  coefficient  matrices 

REM 

REM 

REM 

REM  COEFFICIENT  MATRICES  OF  THE  BASIC  EQUATIONS 

REM 

REM  dX/dT  =  A*x  +  B*U 

REM 

REM  Q  =  C*x  +  D*U 

REM 

REM  where  >  is  the  vector  of  state  variables  and 

REM  u  is  the  vector  of  inputs 

REM 

REM  Matrix  A 

FOR  i=l  to  num_state_var 
FOR  j=l  to  num_state_var 

LET  A(i,j)=G(i,j) *CINV ( i) 

NEXT  j 

FOR  k=l  to  num_state _var  +  num_inputs 
LET  A (i , i)  =  A(i, i)  -  G  (  i  ,  k) *CINV ( i) 

NEXT  k 
NEXT  i 

REM  Matrices  B  and  C 

FOR  i=l  to  num_state_var 
FOR  j=l  to  num_inputs 

LET  B ( i , j ) =G ( i , j  +num_state_var) *CINV ( i ) 


98 


LET  C( j , i)=G(i, j+num_state_var) 

NEXT  j 
NEXT  i 

REM  Matrix  D 

FOR  i=l  to  num_inputs 
FOR  j=l  to  num_inputs 

LET  D ( i , j ) =G ( i+num_state_var , j  +num_state_var ) 
NEXT  j 

FOR  k=l  to  num_state_var  +  num_inputs 

LET  D(i,i)=D(i,i) -G ( i+num_state_var , k) 

NEXT  k 
NEXT  i 

CALL  TF 

CALL  STEADY_STATE_SCLN 
SUB  TF 

REM  Calculation  of  Transfer  Function 


CALL  CALCPHI 
CALL  CALCGAMMAS 
CALL  CALCSCOEFF 
END  SUB 


SUB  CALCPHI 

REM  Calculate  exponential  matrix  PHI 

REM  del=size  of  time  step=l  hour 


99 


LET  del=l . 0 
MAT  ADEL=del*A 

REM  Calculate  the  matrix  row  norm 

FOR  i  =  1  to  num_state_var 
LET  test=0 

FOR  j=l  to  num_state_var 

LET  test  =  test  +  ABS (ADEL(i, j) ) 

NEXT  j 

LET  matrix_row_norm=MAX (matrix_row_norm, test) 

NEXT  i 

REM  Find  the  smallest  integer  such  that  2 "integer  is 
greater  than  or  equal  to 
REM  the  matrix  row  norm 

DO  WHILE  2 *count<matrix_row_norm 
LET  count=count+l 
LOOP 

REM  Divide  matrix  ADEL  by  2 "integer 
MAT  AADJ=(2‘ (-count) ) *ADEL 

REM  Calculate  matrix  row  norm  for  adjusted  ADEL  matrix, 
AADJ 

FOR  i=l  to  num_state_var 
LET  test=0 

FOR  j=l  to  num_state_var 

LET  test  *  test  +  ABS (AADJ (i, j) ) 
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NEXT  j 

LET  aadj_matrix_row_noi.m=MAX  (aadj_matrix_row_norm,  test) 
NEXT  i 

REM  Calculate  number  of  terms  to  keep 

LET  how_far=MIN( ( (3*adj_matrix_row_norm) +6) ,100) 

REM  Calculate  the  exponential  of  the  adjusted  matrix,  AADJ 

MAT  IDEN=IDN(num_state_var) 

LET  num=2 

MAT  PHIADJ=ZER (UBOUND ( A , 1) ,UB0UND(A,2) ) 

DO  WHILE  num<=how_far 
LET  factorial  =1 
FOR  i=num  to  1  STEP  -1 

LET  factorial=factorial*i 
NEXT  i 

LET  factora=( 1/factorial) 

MAT  FACTOR=factora*AADJ 
FOR  j=l  to  num-1 

MAT  FACTOR=FACTOR*AADJ 
NEXT  j 

MAT  PHIADJ=FACTOR+PHIADJ 

LET  num=num+l 

LOOP 

MAT  PHIADJ=AADJ+PHIADJ 
MAT  PHIADJ=IDEN+PHIADJ 

REM  Calculate  exponential  of  matrix  ADEL  from  exponential 
of  matrix  AADJ 
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FOR  i=l  to  count 


MAT  PHINEW=PHIADJ*PHIADJ 
NEXT  i 


IF  count=0  THEN 
MAT  PHI=PHIADJ 
ELSE 

MAT  PHI=PHINEW 
END  IF 

END  SUB 


SUB  CALCGAMMAS 

REM  Calculate  GAMMA1  and  GAMMA2 


MAT  AINV=INV (A) 

LET  determinant=DET 

MAT  TEMP1=PHI-IDEN 
MAT  TEMP1=AINV*TEMP1 
MAT  GAMMA1=TEMP1 *  B 
MAT  TEMP1= ( 1/DEL) *GAMMA1 
MAT  TEMP1=TEMP1-B 
MAT  GAMMA2=AINV*TEMP1 
END  SUB 


SUB  CALCSCOEFF 
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MAT  TEMP1=IDEN*GAMMA2 
MAT  TEMP 1=C* TEMPI 
MAT  S0LD=TEMP1+D 

MAT  WRITE  #3 -.SOLD 


MAT  TABLE=zer  ( num_inputs ,  num_inputs ) 
MAT  TABLE=SCLD 
MAT  RNEW=IDEN 


FOR  i=l  to  num_state_var 
MAT  ROLD=RNEW 
MAT  TEMP1=PHI*R0LD 
LET  trace  =  0 
FOR  j~l  to  num_state_var 

LET  trace=trace+TEMPl ( j , j ) 
NEXT  j 

LET  e--trace/i 

WRITE  #3:  e 

MAT  TEMP1=PHI*R0LD 
MAT  TEMP2=e*IDEN 
MAT  RNEW=TEMP1+TEMP2 
MAT  TEMP1=GAMMA1-GAMMA2 
MAT  TEMP1=R0LD*TEMP1 
MAT  TEMP2=RNEW*GAMMA2 
MAT  TEMP1=TEMP1+TEMP2 
MAT  TEMP1=C*TEMP1 
MAT  TEMP2=e*D 
MAT  SNEW-TEMP1+TEMP2 

MAT  WRITE  #3:  SNEW 
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NEXT  i 


END  SUB 

SUB  PROPERTIES 

REM  Assignment  of  soil  properties 
MAT  CONDUCTIVITY  =  8.640e4 
i J/day-m-K 
MAT  DENSITY  =  1200 
! kg/m* 3 

MAT  SPECIFIC_HT  =  1200  ]J/kg-K 

END  SUB 

SUB  MODEL 

REM  Surface  Description 
REM  Square  Surface 


LET  surface_length=12 

LET  area_of_surface=surface_length‘2  !m*2 

LET  perimeter_of_surf ace=4 *surf ace_length  !m 

LET  c l-area_o f_surf ace/per imeter_of_sur face 

LET  DISTANCE (1) =3. 5  !m 

LET  DISTANCE (2) =14 . 5  I  distance  to  deep  ground  temp 

LET  DISTANCE ( 3 ) =12 . 5  ! distance  to  far  field  temp  m 


LET  DISTANCE (4 )=perimeter_of_surf ace/ (2*pi)  ! distance 

center  to  edge  m 

REM  Model 

REM  Basic  Geometry 


m 
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LET  P(l)=2*pi*( 1.0) /log (DISTANCE (4) /(DISTANCE (4) -1.0) ) 
LET  P(2) =P( 1) 

LET  P(3)=2*pi*4.5*cl 
LET  P ( 4) =P (3) 

LET  H ( 1 ) =DI STANCE ( 1 ) 

LET  H ( 2 ) =DISTANCE ( 2 ) -h ( 1) 

LET  H ( 3 ) =h ( 1 ) 

LET  H ( 4 ) =DISTANCE ( 2 ) -h ( 3 ) 

LET  H(5)=DISTANCE(2)/2 
LET  H (6) =1 . 0 
LET  H (7 ) =0 . 1 

LET  AOV (1) =P (1) *H ( 1 ) 

LET  AOV ( 2 ) =P ( 2 ) *H ( 2 ) 

LET  AOV ( 3 ) =P ( 3 ) *H ( 3 ) 

LET  AOV ( 4 ) =P ( 4 ) *H ( 4 ) 

LET  AOV ( 5 ) =p i * ( (area_of_surface/pi) ‘0. 5-1.0) “2 

LET  AOV (6) =pi* ( (area_of_sur face/pi) ‘0 . 5+2 . 5+0. I5*cl) ~2- 

AOV ( 5 ) 

LET  AOV (7 ) =pi* ( (area_of_sur face/pi) ' 0 . 5+DISTANCE ( 3 ) ) ~2 
AOV (6)  -  AOV ( 5 ) 

LET  VOLUME (1) =AOV(5) *H (5) 

LET  VOLUME (2 )=AOV(6) *H(6) 

LET  VOLUME (3 )=AOV(7) *H(7) 

REM  Set  up  basic  geometry  MATRICES 

MAT  AREA=ZER 

LET  AREA (1,2) , AREA (2,1) =AOV ( 2 ) 

LET  AREA (1,4) , AREA (4,1) =AOV ( 5 ) 

LET  AREA (1,6) , AREA (6,1) = AOV ( 5 ) 
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LET  AREA (2,3) , AREA (3,2) =AOV ( 4 ) 

LET  AREA (2,6) , AREA (6,2) =AOV ( 6 ) 

LET  AREA (2,7) , AREA (7,2) =AOV ( 6 ) 

LET  AREA (3,5) , AREA (5,3) = AO V ( 7 ) 

LET  AREA (3,6) , AREA (6,3) = AO V ( 7 ) 

LET  AREA (4,7) , AREA (7,4) = AOV ( 1 ) 

LET  AREA (5,7) , AREA (7,5) =AOV ( 3 ) 

MAT  LENGTH=ZER 

LET  LENGTH (1,2) , LENGTH (2,1) , LENGTH (4, 7) , LENGTH ( 7 , 4 ) - 
=DISTANCE ( 4 ) 

LET  LENGTH (1,4) , LENGTH (4,1) , LENGTH (2, 7) , LENGTH (7, 2) , LENGTH (3 
, 5 ) , LENGTH (5,3) =DISTANCE ( 1 ) 

LET  LENGTH (1,6) , LENGTH (6,1) , LENGTH (2, 6) , LENGTH (6, 2) , LENGTH (3 
, 6 ) , LENGTH (6,3) =DISTANCE ( 2 ) -DISTANCE ( 1 ) 

LET  LENGTH (2, 3) , LENGTH ( 3 , 2 ) , LENGTH (5, 7) , LENGTH (7 , 5) - 
=DISTANCE ( 3 ) 

END  SUB 


SUB  STEADY_STATE_SOLN 

REM  Calculate  the  steady-state  solution 
DIM  X ( 3 , 1) , CBA (4,4) , GG ( 4 , 4 ) , T ( 4 , 1 ) , Q ( 4 , 1) 
MAT  AINV=INV (A) 

MAT  X=AINV*B 

MAT  CBA=C*X 
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MAT  GG=D-CBA 

MAT  GG=(1/ (3600*24) )*GG 

MAT  WRITE  #3:  GG 

END  SUB 

END 
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APPENDIX  C  TRUBASIC  PROGRAM  QCALC 


TRUE  BASIC  PROGRAM  QCALC 

FOR  USING  MULTIPLE-INPUT  GROUND  TRANSFER  FUNCTIONS 
AND  SCALAR  CONSTANTS  TO  CALCULATE  HEAT  FLUX 


INPUT  PROMPT  "GTF  run  id?  ":gtfrunid$ 

INPUT  PROMPT  "LOCATION  run  id?  " : locid$ 

INPUT  PROMPT  "RUN  NUMBER  id?  ":numid$ 

INPUT  PROMPT  "Number  of  timesteps  to  calculate?" : endtime 


LET  runid$=locid$&numid$ 

LET  outname$=gtfrunid$&runid$&"OUT" 

OPEN  #3: NAME 

"D: \TRUBASIC\INDATA\"&gtfruniu$&"GTFS" , ORGANIZATION 
RECORD, RECSIZE  64 

OPEN  #4: NAME  "D: \TRUBASIC\INDATA\"&locid$&"GSTEMP" , CREATE 
NEWOLD, ORGANIZATION  RECORD, RECSIZE  64 

OPEN  #5 :NAME  "D: \TRUBASIC\DATA\"&outname$ , CREATE  NEWOLD, OR¬ 
GANIZATION  RECORD, RECSIZE  64 
SET  #5: POINTER  BEGIN 
OPEN  #9: NAME 

"D: \TRUBASIC\INDATA\" &locid$&"BTEMP" , ORGANIZATION 
RECORD, RECSIZE  64 
RESET  #9: BEGIN 

READ  #3 :gtfrundate$ , gtfruntime$ , gtfrunid$ 

DIM  SO (4 , 4 ) , SI (4 , 4 ) , S2 ( 4 , 4 ) , S3 ( 4 , 4 ) ,E(0  TO  3) 

DIM  SSTF (4,4) 

DIM  QSTART ( 4 ) , TSTART ( 4 ) 

DIM  Q (4 , 0  TO  3 ) , T ( 4 , 0  TO  3) 
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DIM  QO ( 4 ) , TO ( 4 ) , QSAVE ( 4 ) 

DIM  FACTORO (4,0  TO  3) , FACT0R1 (4,0  TO  3 ) , FACTOR2 (4,0  TO 
3 ) , FACTOR3 (4,0  TO  3),FACTORE(0  TO  4) 

DIM  QDIFF (  4 )  , QCONV (4,0  TO  3) , PDIFF(4) ,QDIFFCONV(4) 

WRITE  #5:  outname$,gtfrunid$ 

WRITE  #5: end time 

REM  Initialize  variables 

LET  E (0) =0 . 0 
LET  num_inputs=4 
LET  num_state_var=3 
LET  iteration  =0 
LET  er=0 . 001 
LET  timestep=0 

REM  Read  in  Multiple-input  GTF  coeffcients  and  scalar  con¬ 
stants 

MAT  READ  #3: SO 
READ  #3:  E ( 1) 

MAT  READ  #3:  SI 
READ  #3:  E (2) 

MAT  READ  #3:  S2 
READ  #3:  E (3 ) 

MAT  READ  #3:  S3 
MAT  READ  #3:  SSTF 
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REM  Set  up  starting  temperture  and  flux  matrices 


RESET  #4:  RECORD  2 
READ  #4 : TSTART ( 2 ) 

READ  #9: TSTART (1) 

RESET  #9: BEGIN 

INPUT  PROMPT  "DEEP  GROUND  TEMP  [K]:  " : TSTART ( 3 ) 

LET  TSTART (4 )=TST ART (1) 

RESET  #4: RECORD  2 
MAT  QSTART=SSTF*TSTART 
MAT  QSAVE=ZER 
FOR  inp=l  to  num_inputs 
FOR  hr=0  to  num_state_var 
LET  Q ( inp, hr) =QSTART ( inp) 

LET  T ( inp, hr) =TSTART ( inp) 

NEXT  hr 
NEXT  inp 

MAT  Q= (3600*24) *Q 

DO 

REM  If  system  is  not  yet  initialized,  initialize  system 


IF  timestep=0  THEN 
CALL  INITLOOP 
LET  timestep=timestep+l 
END  IF 

IF  timestep=l  THEN 
RESET  #4: RECORD  2 
END  IF 
CALL  CAL CQ 

LET  timestep=timestep+i 
LOOP  UNTIL  timestep=>endtime 
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SUB  INITLOOP 

REM  Initialize  system  by  iterating  on  first  30  days  of  data 
DO 

RESET  #4: RECORD  2 

LET  iteration=iteration+l 

FOR  in=l  to  num_inputs 

LET  QDIFF ( in) =Q ( in, 0) -QSAVE ( in) 

IF  QSAVE ( in) <>0 . 0  THEN 

LET  PDIFF ( in) » (QDIFF ( in) / (ABS (QSAVE (in) ) ) ) *100 
END  IF 

LET  QSAVE(in)=Q(in, 0) 

NEXT  in 
LET  cycle=0 
DO  WHILE  cycle<=30 
LET  cycle=cycle+l 

MAT  QDIFFCONV=(l/ (3600*24) ) *QDIFF 
FOR  inp  =1  to  num_inputs 

FOR  hour=num_state_var-l  to  0  step  -1 
LET  Q ( inp, hour+l)=Q( inp, hour) 

LET  T ( inp, hour+l)=T( inp, hour) 

NEXT  hour 
NEXT  inp 
IF  END  #4  THEN 

RESET  #4: record  2 
END  IF 

IF  END  #9  THEN 
RESET  #9: BEGIN 
END  IF 

READ  #4:  T (2 , 0) 

LET  T ( 3 , 0 ) =TSTART ( 3 ) 
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REM  Set  slab  edge  temperature  =  slab  center  temperature 


LET  TSTART ( 4 ) =TSTART ( 1 ) 

CALL  CALCQ 
LET  TO ( 1 ) =T (1,0) 

LET  TO (2 ) =T (2,0) 

LET  TO ( 3 ) =T (3,0) 

LET  TO (4) =T (4,0) 

LET  QO (1) =Q (1,0) 

LET  QO ( 2 ) =Q (2,0) 

LET  QO ( 3 ) =Q (3,0) 

LET  QO ( 4 ) =Q (4,0) 

MAT  Q0=(1/ (3600*24) ) *Q0 
LOOP 

REM  If  system  does  not  initialize  in  200  iterations, 
nate  initialization 

IF  iteration  >  200  THEN 
EXIT  DO 
END  IF 

LOOP  UNTIL  ABS ( PDIFF ( 1 ) ) <er  AND  ABS ( PDIFF ( 2 ) ) <er  AND 
ABS (PDIFF (3) ) <er  AND  ABS ( PDIFF ( 4 ) ) <er  AND  iterations 
END  SUB 

SUB  CALCQ 

IF  timestep>0  THEN 

FOR  inp  =1  to  num_inputs 

FOR  hour=num_state_var- J  to  0  step  -1 
LET  Q ( inp, hour+1 ) =Q ( inp, hour) 

LET  T(inp,hour+l)=T(inp,hour) 

NEXT  hour 


termi- 


112 


LET  Q(inp,0)=0.0 
NEXT  inp 

IF  END  #4  THEN 

RESET  #4: record  2 
END  IF 

IF  END  #9  THEN 

RESET  #9: BEGIN 
END  IF 

READ  #4:  T (2 , 0) 

LET  T  (  3 , 0 ) =TSTART ( 3 ) 

READ  #9:  T(1,0) 

LET  TSTART ( 4 ) =TSTART ( 1 ) 

END  IF 

MAT  FACTORO=SO*T 
MAT  FACTORIES 1*T 
MAT  FACTOR2=S2*T 
MAT  FACTOR3=S3*T 

LET  factor_format$  =  "-#.########## - " 

LET  qsum=0.0 

FOR  L=1  TO  4 

LET  FACTORE ( L) =E ( 1 ) *Q ( L, 1 ) +E ( 2 ) *Q (L, 2 ) +E ( 3 ) *Q (L, 3 ) 

LET  Q (L, 0) =FACTORO (L, 0) +FACTOR1 (L, l)+FACTOR2 ( L, 2 ) +FACTOR3 
(L, 3) -FACTORE (L) 

LET  qsum=qsum+Q (L, 0) 

NEXT  L 

MAT  QCONV  =  (1/(3600*24) ) *Q 
LET  qsumconv= (1/ (3600*24) ) *QSUM 
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IF  timestep=0  THEN 

LET  q_format$="“# .########  « 

ELSE 

WRITE  #5 :QCONV (1,0) ,QCONV(2,0) ,QCONV(3,0) ,QCONV(4,0) 
END  IF 

END  SUB 

END 
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